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Abstract 

In  recent  years,  imaging  through  atmospheric  turbulence  has  interested  mil¬ 
itary  scientists  seeking  to  improve  optical  surveillance  of  satellites  for  intelligence 
gathering  purposes.  Adaptive  optics  was  a  step  toward  achieving  diffraction-limited 
resolution  from  ground-based  telescopes.  Unfortunately,  adaptive  optics  only  par¬ 
tially  compensate  for  atmospheric  blurring,  therefore  post-processing  of  images  is 
required.  Processing  methods  in  use  today  require  knowledge  of  the  impulse  re¬ 
sponse  of  the  optical  system  to  reconstruct  imagery,  but  this  information  is  seldom 
known.  This  thesis  looks  at  a  new  method  of  processing  compensated  imagery,  called 
blind  deconvolution,  which  assumes  very  little  or  no  a  •priori  information  about  the 
impulse  response.  In  particular,  this  investigation  analyzes  the  performance  of  Lane’s 
unconstrained  mininoization  method  of  blind  deconvolution  as  modified  by  Jefferies 
and  Christou.  The  unconstrained  minimization  technique  is  applied  to  simulated 
single  and  binary  star  images  corrupted  by  photon  noise.  Results  reveal  that  prior 
knowledge  of  the  cutoff  frequency  of  the  system  greatly  enhances  the  ability  of  the 
algorithm  to  achieve  accurate  estimates  of  the  object  when  measurements  contain 
relatively  few  photo  events.  Additionally,  this  study  discovered  that  estimates  are 
highly  dependent  upon  the  choice  of  the  support  region.  Analysis  also  shows  that 
the  algorithm  produces  estimates  containing  frequency  content  above  the  diffraction- 
limit.  The  presence  of  this  high-frequency  information  may  invalidate  this  method 
as  a  useful  means  to  reconstruct  imagery. 
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EFFECTS  OF  PHOTON  NOISE  ON  UNCONSTRAINED 
MINIMIZATION  TECHNIQUES  FOR  ITERATIVE 
BLIND  DECONVOLUTION 


/.  Introduction 

1.1  Motivation 

For  centuries,  scientists  and  philosophers  alike  have  looked  to  the  stars  for 
numerous  purposes:  for  meaning  in  life  as  well  as  a  better  understanding  of  our 
universe  and  our  planet.  The  invention  of  the  telescope  permitted  man  a  closer  look 
at  the  stars,  but  many  soon  realized  that  the  starHght  they  beheld  was  distorted  by 
the  atmosphere.  Isaac  Newton  surmised  in  1704  that  better  observations  might  be 
made  on  the  peaks  of  high  mountains  to  mitigate  “the  confusion  of  the  rays  which 
arises  from  tremors  of  the  atmosphere”  (16). 

Today,  the  telescope  has  uses  far  beyond  what  Newton  and  his  contemporaries 
dreamed.  The  launch  of  Sputnik  in  1957  threw  the  United  States  into  a  flurry  of 
activity  to  respond  in  kind  while  apprehension  over  the  Soviet’s  ability  to  place  an 
object  in  orbit  created  an  entirely  new  endeavor — space  surveillance.  The  Air  Force 
met  this  new  mission  requirement  by  first  employing  the  Baker-Nunn  satellite  track¬ 
ing  camera  at  five  locations  around  the  world  (15).  The  Baker-Nunn  camera  used  a 
modified  Schmidt  telescope  with  excellent  resolution  to  photograph  and  identify  ob¬ 
jects  in  space.  Since  the  early  1960’s,  the  space  surveillance  mission  has  shifted  from 
visual  optics  to  radar  which  is  not  limited  to  operations  at  night.  However,  optical 
surveillance  maintains  a  key  role  in  accomphshing  the  space  surveillance  mission. 
The  ground-based  electro-optical  deep  space  surveillance  system  (GEODSS)  sensors 
are  responsible  for  collecting  over  65  percent  of  the  tracking  data  for  deep  space  ob- 
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jects  (15).  According  to  the  Air  Force’s  Space  Handbook,  space  surveillance  includes 
“the  ability  to  surveil  and  monitor  continuously  all  significant  mihtary  activities  in 
space”  (15).  Today,  Air  Force  leaders  not  only  want  to  know  the  satellite  and  its 
owner,  but  also  what  missions  it  performs  and  what  equipment  it  carries.  This  re¬ 
quires  much  more  resolution  than  current  GEODSS  telescopes  or  even  phased-array 
radar  can  provide. 

Research  toward  improving  this  resolution  is  accomplished  at  the  U.  S.  Air 
Force  Maui  Optical  Station  (AMOS).  AMOS  houses  one  of  the  most  advanced  tele¬ 
scopes  in  the  world  for  the  purpose  of  imaging  objects  in  space.  It’s  1.6  meter 
telescope  obtains  short-exposure,  high  resolution  imagery  through  the  use  of  an 
adaptive  optics  system.  Since  adaptive  optics  only  partially  compensates  for  the 
turbulence-induced  aberrations,  AMOS  also  rebes  on  post-processing  techniques  to 
improve  overall  image  quality  (23). 

Much  work  has  gone  into  developing  algorithms  and  techniques  which  can  im¬ 
prove  the  quality  of  blurred  images.  Most  techniques  require  an  understanding  of  the 
impulse  response,  also  called  the  point  spread  function  (PSF),  and  some  knowledge 
of  the  statistical  nature  of  the  noise  inherent  in  the  optical  system.  Adaptive  optics 
images  result  from  the  actual  object  and  the  compensated  PSF,  which  includes  ef¬ 
fects  from  atmospheric  turbulence  not  corrected  by  the  adaptive  optics  system.  To 
alleviate  the  problem  of  the  unknown  PSF,  one  image  processing  method,  blind  de- 
convolution,  makes  successive  guesses  at  both  the  object  and  point  spread  function 
such  that  they  converge  on  the  “true”  object  and  PSF. 

Blind  deconvolution  provides  the  ability  to  make  estimates  of  both  an  object 
and  a  PSF  without  a  priori  knowledge  of  the  blurring  function.  Though  several 
papers  have  been  written  on  the  topic  of  blind  deconvolution,  the  work  completed 
by  R.  G.  Lane  (13)  can  most  readily  be  applied  to  short  exposure  compensated 
imagery.  The  Lane  method  as  modified  by  S.  M.  Jefferies  and  J.  C.  Christou  allows 
for  great  flexibility  since  it  makes  estimates  of  both  the  object  and  the  PSF  given 
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any  amount  of  information  available  on  the  impulse  response — from  no  information 
to  complete  knowledge  of  the  optical  system  (12). 

Since  the  adaptive  optics  system  at  AMOS  takes  short-exposure  images  con¬ 
taining  very  httle  light,  any  post-processing  technique  must  have  the  ability  to  recon¬ 
struct  images  in  the  presence  of  photon  noise.  Jefferies  and  Christou  only  addressed 
the  issue  of  additive  noise.  Thus,  before  recommending  this  technique  to  the  opera¬ 
tional  users  on  Maui,  the  question  should  be  asked:  How  well  does  the  unconstrained 
minimization  algorithm  perform  for  photon  hmited  images? 

1.2  Problem  Statement 

This  thesis  seeks  to  document  the  abihty  of  the  iterative  deconvolution  al¬ 
gorithm  to  reconstruct  images  in  the  presence  of  photon  noise.  The  algorithm’s 
performance  is  evaluated  as  a  function  of  differing  hght  levels  and  differing  amounts 
of  prior  information  concerning  the  blurring  function  and  the  noise. 

1.3  Approach 

To  accompUsh  the  task  of  evaluating  the  unconstrained  minimization  tech¬ 
nique,  several  things  were  put  in  order.  First,  S.  M.  Jefferies  and  J.  C.  Christou 
gave  permission  to  use  their  FORTRAN  code  which  implements  the  unconstrained 
minimization  technique  through  their  iterative  deconvolution  algorithm  (IDA)  as 
discussed  in  their  1993  paper  (12).  Second,  a  model  had  to  be  developed  which 
would  isolate  photon  noise  from  all  other  noise  and  blurring  effects.  Finally,  experi¬ 
mentation  was  designed  around  answering  the  following  questions. 

1.  How  can  photon  noise  be  modeled  for  computer  simulation? 

2.  At  what  Hght  level  can  IDA  no  longer  improve  the  measured  image? 

3.  Do  the  modifications  made  by  Jefferies  and  Christou  provide  significant 
improvement  over  the  Lane  method  for  reconstructing  photon  Hmited  imagery? 
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4.  How  can  prior  information  regarding  the  point  spread  function’s  cutoff 
frequency  best  be  utilized  in  the  reconstruction? 

As  a  result,  computer  simulation  was  selected  as  the  best  means  of  both  cre¬ 
ating  photon  hmited  images  and  testing  the  bUnd  deconvolution  algorithm. 

1.4  Scope 

This  thesis  seeks  to  simulate  post-processing  of  images  of  exo-atmospheric  ob¬ 
jects  obtained  through  an  adaptive  optics  system  and  specifically  analyze  the  effects 
of  photon  noise  in  the  measured  images.  The  single  simulated  images  are  Hmited  to 
point  sources  of  differing  intensity  and  separation.  These  images  are  directly  analo¬ 
gous  to  single  and  binary  star  images.  Multiple  or  composite  images  and  extended 
objects,  such  as  sateUites,  are  not  addressed  within  this  study. 

Jefferies  and  Christou  modified  the  Lane  method  by  adding  band-Hmit  and 
Fourier  modulus  error  terms  to  the  objective  function  for  minimization.  Since  the 
IDA  code  allows  selective  use  of  both  of  these  error  metrics,  only  the  band-Hmit 
metric  was  tested.  The  following  analysis  does  not  include  use  of  the  Fourier  modulus 
error  metric,  since  its  appHcation  requires  Fourier  modulus  information  obtained 
through  speckle  interferometry  (13).  One  goal  of  this  paper  is  to  reveal  the  utiHty 
of  pure  bHnd  deconvolution,  thus  as  Httle  additional  information  as  possible  is  used 
in  the  reconstruction  process. 

Additionally,  this  thesis  Hmits  its  scope  to  post-processing  of  detected  images 
only.  It  does  not  discuss  design  of  an  adaptive  optics  system.  Optical  design  issues 
such  as  properties  of  optical  materials  or  detector  array  functions  are  not  addressed 
here.  FinaUy,  this  paper  does  not  consider  real-time  processing  of  images  nor  does 
it  specifically  cite  the  computational  hardware  necessary  to  process  images. 
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1.5  Summary  of  Key  Results 

Several  results  were  obtained  in  this  analysis  which  weigh  heavily  on  the  fu¬ 
ture  study  of  this  algorithm  to  process  images  of  extended  objects.  Proof  that  the 
unconstrained  minimization  technique  will  work  for  extended  objects  is  absolutely 
necessary  before  any  proposal  can  be  made  regarding  its  operational  use  at  AMOS. 
The  key  results  obtained  by  this  study  are  as  follows: 

1.  Support  size  and  shape  greatly  affect  the  reconstruction.  Incorporating 
knowledge  of  the  shape  and  especially  the  location  of  the  first  ring  of  the  actual 
point  spread  function  in  the  support  arrays  produces  object  estimates  much  closer 
to  the  “true”  object. 

2.  Information  regarding  the  diffraction-hmited  spatial  frequency  cutoff,  or 
band-limit,  of  the  measured  image  is  necessary  to  reconstruct  any  photon  limited 
image,  even  when  the  amount  of  Hght  present  is  relatively  high.  Without  the  band- 
limit  input,  the  noise  in  the  measured  image  is  multiplied  in  the  estimates  produced 
by  IDA  similar  to  noise  effects  in  an  inverse  filter. 

3.  The  IDA  code  performs  well  for  very  high  hght  levels.  As  the  hght  level  is 
reduced,  Poisson  noise  effects  due  to  the  random  arrival  of  photo-events  become  so 
strong  that  single  images  at  low  hght  levels  are  almost  impossible  to  restore  correctly. 

4.  Unconstrained  minimization  techniques  allow  for  the  creation  of  information 
above  the  cutoff  frequency — an  attempt  at  super-resolution.  The  high  frequency 
content  in  the  estimated  object  spectrum  and  optical  transfer  function  (OTF)  leads 
one  to  question  the  vahdity  of  the  iterative  deconvolution  algorithm,  especially  since 
the  point  spread  function  is  band-hmited  due  to  the  physical  nature  of  the  optics. 

1.6  Chapter  Outlines 

A  brief  synopsis  of  the  remainder  of  the  thesis  foUows. 
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1.6.1  Chapter  Two.  Chapter  2  reviews  some  earlier  work  in  the  area  of 
blind  deconvolution  and  outlines  the  algorithm  derived  by  Lane  and  modified  by 
Jefferies  and  Christou. 

1.6.2  Chapter  Three.  Chapter  3  lays  down  the  terminology,  the  compari¬ 
son  metrics  and  the  objectives  which  allow  for  computer  simulation  of  the  iterative 
deconvolution  algorithm. 

1.6.3  Chapter  Four.  Chapter  4  details  three  experiments  designed  to  test 
the  blind  deconvolution  technique  and  analyzes  the  results  of  the  experimentation 
with  the  main  focus  on  the  effects  of  photon  noise. 

1.6.4  Chapter  Five.  Chapter  5  summarizes  the  important  results  uncovered 
throughout  this  thesis  effort  and  concludes  with  recommendations  for  further  study 
on  extended  objects. 

1.6.5  Appendix  A.  Appendix  A  documents  the  iterative  deconvolution 
algorithm  (IDA)  code  as  modified  for  this  thesis.  This  appendix  provides  consider¬ 
able  detail  on  the  use  of  the  FORTRAN  code  provided  by  S.  M.  Jefferies  and  J.  C. 
Christou  which  performs  blind  deconvolution  through  unconstrained  minimization. 
Included  in  the  appendix  are  the  appropriate  inputs  used  by  the  program  and  the 
proper  format  for  the  different  data  files.  This  appendix  also  contains  some  helpful 
hints  regarding  modification  of  the  code  for  future  study. 

1.6.6  Appendix  B.  Appendix  B  contains  a  compilation  of  the  simulated 
input  images  utilized  in  three  experiments. 

1.6.7  Appendix  C.  Appendix  C  consists  of  estimates  produced  by  the  Lane 
method  in  Experiment  1  for  noise-free  images. 
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1.6.8  Appendix  D.  Appendix  D  compiles  the  object  estimates  made  during 
Experiment  2  for  photon  hmited  images  using  both  the  strict  Lane  method  and  the 
modified  algorithm  of  Jefferies  and  Christou. 

1.6.9  Appendix  E.  Appendix  E  fists  the  output  estimates  found  in  Exper¬ 
iment  3  utilizing  the  band-limit  error  term  in  unconstrained  minimization. 
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II.  Background 


2. 1  Introduction 

This  chapter  presents  the  basic  framework  necessary  to  understand  the  funda¬ 
mental  problem  of  atmospheric  turbulence,  its  relation  to  adaptive  optics,  and  why 
post-processing  of  images  is  necessary  after  adaptive  optics  correction.  In  addition 
to  the  fundamentals  of  atmospheric  imaging,  background  material  is  presented  on  a 
particular  post-processing  technique  called  bhnd  deconvolution.  After  a  brief  survey 
of  various  bhnd  deconvolution  algorithms,  this  chapter  concludes  with  a  detailed 
look  at  the  unconstrained  minimization  method  derived  by  R.  G.  Lane  and  modified 
by  S.  M.  Jefferies  and  J.  C.  Christou. 

2.2  Atmospheric  Turbulence 

The  resolution  of  images  formed  by  large  optical  telescopes  is  directly  at¬ 
tributable  to  atmospheric  turbulence.  According  to  the  “turbulent  eddy”  model, 
this  degradation  is  due  to  random  inhomogeneities  in  the  refractive  index  of  air  (8). 
Differential  heating  of  the  Earth’s  surface  creates  this  phenomenon,  producing  large 
scale  temperature  gradients.  Convection  and  turbulent  wind  flow  break  up  these 
large  scale  inhomogeneities  into  smaller  scale  “eddies”.  Each  turbulent  eddy  has 
a  unique  refractive  index,  which  modulates  the  amphtude  and  phase  of  a  propa¬ 
gating  wavefront  both  temporally  and  spatially.  Amphtude  modulation  results  in 
scintillation  as  observed  in  the  twinkhng  of  stars,  whereas  phase  modulation  results 
in  random  image  motion  (known  as  tilt)  and  phase  aberration.  For  the  case  of 
ground-based  imaging  of  exo-atmospheric  objects,  the  effects  of  phase  modulation 
are  generally  more  severe  than  amphtude  modulation  (5). 

Astronomers  concerned  with  the  propagation  of  hght  waves  through  turbulence 
have  developed  various  parameters  for  characterizing  the  severity  of  image  degrada¬ 
tion  due  to  turbulence.  Astronomers  typically  use  such  parameters  to  compare  the 
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relative  seeing  quality  of  candidate  sites  for  new  observatories.  One  of  the  most 
convenient  and  widely  used  measures  of  seeing  quality  was  introduced  by  Fried  (6) 
and  is  denoted  tq.  It  is  defined  as  the  effective  diameter  of  a  telescope  for  which  the 
Strehl  resolution  of  the  telescope  is  equal  to  the  Strehl  resolution  associated  with 
the  atmospheric  optical  transfer  function  (OTF)  (18).  ro  is  a  function  of  the  zenith 
angle  of  the  path  of  propagation,  the  wavelength,  and  the  turbulence  strength  (18). 
Since  the  effects  of  turbulence  on  optical  propagation  are  random  in  nature,  ro  is 
also  random.  Typical  values  for  ro  at  a  good  observatory  range  from  5  cm  for  mod¬ 
erately  poor  seeing  to  20  cm  for  exceptional  seeing  (8).  In  addition  to  its  practical 
use  as  a  measure  of  relative  seeing  quality,  ro  is  widely  used  in  expressions  for  the 
atmospheric  OTF  (8).  Since  it  acts  as  an  effective  telescope  diameter,  ro  simphfies 
the  form  of  the  OTF,  and  aids  in  understanding  the  effect  of  the  atmosphere  relative 
to  the  true  diameter  of  the  optics. 

Since  the  atmospheric  turbulence  changes  randomly  as  a  function  of  both  time 
and  space,  no  two  images  taken  through  the  atmosphere  can  have  the  same  atmo¬ 
spheric  OTF  unless  they  are  imaged  through  the  same  part  of  the  atmosphere  at  the 
exact  same  time.  The  usual  method  of  modeling  the  imaging  process  is  by  assuming 
the  process  is  hnear  and  shift  invariant  (2).  Given  an  object,  f{x,y),  and  a  point 
spread  function  (PSF),  h{x,y),  where  {x,y)  is  a  position  vector  locating  an  arbitrary 
point  in  image  space,  the  measured  image,  g{x,y),  is  given  by 

9{^,y)  =  (2-i) 

and  by  the  convolution  theorem, 

g{x,y)  ^  G{u,v)  =  F{u,v)H{u,v),  (2.2) 

where  denotes  the  two-dimensional  convolution  of  the  two  functions,  the  capital 
letters  represent  the  Fourier  transforms  of  the  respective  lowercase  functions  and  the 
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identifies  a  Fourier  transform  pair.  Points  in  Fourier  space  are  given  by  the 
position  vector  {u,v).  Therefore, 

/+00  ^+00 

/  f{x,y)  exp[j2TT{ux  +  vy)]dxdy,  (2.3) 

-00  J —  00 

where  F[u,v)  represents  the  Fourier  transform  of  f[x,y). 

If  the  object  can  be  assumed  to  be  stationary,  then  one  only  needs  to  determine 
under  what  conditions  the  PSF,  h(x,y),  is  shift  invariant  to  have  a  linear,  shift 
invariant  system.  Since  the  PSF  is  the  inverse  Fourier  transform  of  the  optical 
transfer  function,  H(u,v),  the  PSF  varies  both  spatially  and  temporally  as  well.  To 
freeze  the  time- varying  shift  in  the  PSF  or  OTF,  short-exposure  images  on  the  order 
of  miUi-seconds  are  required.  To  ensure  shift  invariance,  each  object  must  be  viewed 
through  a  very  small  cone  of  the  atmosphere.  This  small  angle  is  known  as  the 
isoplanatic  angle.  Equation  2.1  above  is  valid  only  for  short-exposures  and  when  the 
object  lies  within  the  isoplanatic  angle. 

2.3  Adaptive  Optics 

Though  atmospheric  imaging  requires  short-exposures  to  be  imaged  through 
very  small  solid  angles,  real-time  compensation  of  the  random  phase  modulation  due 
to  turbulence  can  be  achieved  through  the  use  of  adaptive  optics.  To  accomplish 
this,  the  adaptive  optics  system  must  measure  these  phase  aberrations  and  compen¬ 
sate  for  them  almost  instantaneously  using  a  deformable  mirror.  Determining  the 
phase  perturbation  requires  a  wave-front  sensor  which  observes  a  near-by  natural 
star  or  a  synthetic  laser  guide  star  (21).  The  wave-front  sensor  communicates  the 
phase  aberrations  present  in  a  particular  isoplanatic  region  of  the  atmosphere  at  a 
specific  time  to  the  electronics  which  deform  the  mirror  appropriately.  In  the  AMOS 
system,  a  monolithic  piezoelectric  mirror  (MPM)  has  168  actuators  which  push  or 
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pull  on  different  locations  around  the  mirror  to  correct  for  phase  errors  caused  by 
atmospheric  turbulence  (19). 


2.^  Post-processing  Techniques 

Unfortunately,  residual  phase  errors  still  remain  even  after  adaptive  optics 
correction  (21).  Therefore,  post-processing  of  images  is  required  to  correct  for  the 
effects  of  residual  phase  errors  and  to  sharpen  the  image  since  high  spatial  frequency 
information  is  attenuated  by  the  OTF.  Sharpening  of  the  image  can  be  accomphshed 
through  deconvolution  of  the  point  spread  function  from  the  measured  image.  Quite 
often  an  ensemble  of  images  are  averaged  and  the  average  point  spread  function  can 
be  deconvolved  from  the  average  measured  image  (19).  Roggemann  [1992]  uses  both 
a  modified  inverse  filter  and  a  pseudo-Wiener  filter  to  make  estimates  at  an  object 
spectrum  given  an  ensemble  of  image  measurements  from  an  adaptive  optics  system 
(19). 

Both  methods  discussed  by  Roggemann  require  some  knowledge  of  the  OTF 
and  the  noise  present  in  the  images.  If  {G{u,v))  represents  the  average  spectrum 
of  the  measured  images,  then  the  object  spectrum  can  be  estimated  by  the  inverse 


F{u,v) 


{G{u,v)) 

H{u,v) 


and  the  pseudo-Wiener  filter, 


T(u,v)  _  r  a  i‘ 


where  H{u,v)  is  an  estimate  of  the  average  compensated  OTF,  SNR(ii,u)  is  an 
estimate  of  the  image  spectrum  signal-to-noise  ratio  and  f3  is  &  parameter  of  choice 
(19).  Errors  in  the  object  spectrum  estimate,  F{u,v),  occur  any  time  H{u,v),  or 
SNR(u,  u)  are  not  chosen  exactly.  Since  H(u,v)  in  particular  is  seldom  known,  new 


2-4 


techniques  have  recently  been  proposed  which  estimate  both  the  object  spectrum 
and  the  OTF. 

2.5  Survey  of  Blind  Deconvolution  Techniques 

T.  G.  Stockham  coined  the  term  “blind  deconvolution”  as  signal  or  image 
recovery  from  a  large  collection  of  differently  blurred  versions  of  the  same  signal  or 
image  (22).  It  is  interesting  to  note  that  one  of  Stockham’s  apphcations  for  bhnd 
deconvolution  involved  the  hopeM  recovery  of  pristine  vocalizations  of  famous  opera 
singers  from  old  gramophone  records  (2).  Today,  blind  deconvolution  of  imagery 
categorizes  the  techniques  which  seek  to  estimate  an  object’s  intensity  directly,  either 
from  a  single  image  or  ensemble  of  images  without  prior  knowledge  of  the  PSF.  Most 
techniques  involve  an  iterative  process  which  starts  with  an  initial  estimate  at  the 
object  and  PSF  and  converges  toward  the  “true”  object  and  PSF  with  each  iteration. 
These  methods  require  a  vast  amount  of  computer  processing  needing  hundreds  to 
thousands  of  iterations  to  complete. 

The  blind  deconvolution  problem  has  been  approached  several  different  ways. 
The  first  method  involves  an  iterative  process  in  which  the  object  and  PSF  estimates 
are  successively  modified  by  certain  image  plane  constraints.  A  second  approach 
involves  maximum-hkehhood  estimation  using  the  expectation-maximization  algo¬ 
rithm.  The  third  method  formulates  the  problem  as  an  unconstrained  minimization 
problem  and  uses  a  conjugate  gradient  minimization  algorithm  to  iteratively  deter¬ 
mine  the  object  and  PSF.  The  proceeding  sections  briefly  sketch  each  of  the  three 
approaches  to  the  blind  deconvolution  problem. 

2.5.1  Ayers-Dainty  method.  In  1988,  G.  R.  Ayers  and  J.  C.  Dainty  pub¬ 
lished  the  first  paper  containing  a  working  blind  deconvolution  algorithm  (1).  The 
Ayers  and  Dainty  technique,  shown  pictoriaUy  in  Figure  2.1,  requires  an  initial  guess 
at  the  PSF  and  then  utihzes  an  iterative  method  which  subsequently  constrains  each 
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Figure  2.1  Ayers-Dainty  blind  deconvolution  algorithm.  Positivity  is  applied  at 
steps  4  and  8  while  the  Fourier  domain  constraints  are  appHed  at  steps 
2  and  6. 
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object  and  PSF  estimate  in  the  image  domain  and  the  frequency  domain.  After  an 
initial  guess  at  the  object,  /o(»,y),  is  provided,  the  Ayers-Dainty  algorithm  performs 
the  following  eight  steps  iteratively. 

1.  Fourier  transform  the  object  estimate,  f{x,y),  to  obtain  an  object  spectrum 
estimate,  F(u,v). 

2.  Form  a  new  estimate  of  the  OTF,  H(u,v),  from  the  measured  image  spec¬ 
trum,  G{u,v),  and  F(u,v).  Rather  than  apply  a  simple  inverse  filter,  the  Ayers- 
Dainty  approach  applies  a  Fourier  domain  constraint.  This  constraint  ensures  the 
product  of  both  F  and  H  equal  G  while  preventing  small  errors  in  G  from  resulting 
in  large  errors  in  H  where  F  is  extremely  small.  Also,  the  Fourier  domain  constraint 
forms  a  weighted  average  of  the  previous  estimate  of  the  OTF  and  the  inverse  fil¬ 
tered  estimate,  G/F.  The  Fourier  constraint  applied  in  the  Ayers-Dainty  algorithm 
is  summarized  below. 

If  \G{u,v)\  <  noise  level  , 

Hi+i(u,v)  =  Hi{u,vy,  (2.6) 

if  \F{u,v)\  >  |G(u,i;)|, 

•  =  +  (2.7) 

F[u,v) 

if  \F{u,v)\  <  |G(«,n)|, 

Hi+i{u,v)  Hi{u,v)  ^G{u,vy 

where  0  </?<  1  and  /3  is  set  before  algorithm  is  run  and  i  and  i  -|-  1  represent  the 
current  iteration  and  next  iteration,  respectively  . 

3.  Inverse  Fourier  transform  the  OTF  estimate,  H{u,v),  to  obtain  a  PSF 
estimate,  h{x,y). 
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4.  Impose  image  plane  constraints  on  h[x,y)  to  obtain  h{x,y).  This  image 
plane  constraint  forces  every  pixel  within  the  estimate  to  be  positive  by  setting  each 
negative  pixel  to  zero.  The  algorithm  checks  every  pixel  value  at  each  iteration.  The 
algorithm  keeps  track  of  the  energy,  or  negative  values,  removed  by  summing  their 
absolute  values.  Then  the  algorithm  equally  redistributes  the  removed  energy  to  all 
pixels  in  the  PSF  estimate. 

5.  Fourier  transform  the  constrained  PSF  estimate,  h{x,y),  to  obtain  an  OTF 
estimate,  H(u,v). 

6.  Impose  the  same  Fourier  domain  constraints  Hsted  in  step  2  to  determine  a 
new  value  for  the  object  spectrum,  F{u,v),  by  replacing  H  with  F  and  visa  versa. 

7.  Inverse  Fourier  transform  the  object  spectrum  estimate,  F(u,v),  to  obtain 
an  object  estimate,  f(x,y). 

8.  Impose  the  positivity  constraint  on  f{x,y)  to  obtain  a  new  object  estimate, 
f{x,y).  If  the  new  object  estimate  is  not  adequate,  then  return  to  step  1. 

The  Ayers-Dainty  approach  works  well  for  images  without  noise  corruption. 
However,  the  method  requires  further  refinement  in  order  to  accommodate  contam¬ 
inated  or  complex  images.  In  1989,  B.  L.  K.  Davey,  R.  G.  Lane,  and  R.  H.  T.  Bates 
modified  the  Ayers-Dainty  technique  to  handle  small  amounts  of  noise  in  the  image 
by  implementing  a  Wiener  filter  in  the  Fourier  domain  constraint  to  obtain  each 
consecutive  estimate  of  the  object  spectrum  and  the  optical  transfer  function.  In 
the  image  domain,  Davey,  et  al.  apply  a  support  constraint  in  addition  to  positivity. 
The  support  constraint  restricts  the  object  and  PSF  estimate  pixels  to  zero  outside 
a  specified  region.  This  method  requires  the  input  of  a  binary  mask  which  gives 
the  extent  of  the  object  and  another  mask  containing  known  extent  of  the  PSF  (4). 
Davey’s  modification  of  Ayers-Dainty  revealed  that  blind  deconvolution  was  possible 
for  images  corrupted  by  random  additive  noise. 


2-8 


Another  modification  to  the  Davey  approach  by  Miura  handles  multiple  speckle 
images  of  the  same  object  (14).  Each  individual  image,  gn{x,y)  can  be  modeled  in 
the  following  way 

9r.{x,  y)  =  /(x,  y)  *  K{x,  y),  (2.9) 

where  hn{x,y)  represents  the  nth  instantaneous  PSF  and  the  asterisk,  denotes 
two-dimensional  convolution.  The  Miura  algorithm  performs  bHnd  deconvolution  on 
each  speckle  image  and  averages  the  resulting  object  estimates  after  each  iteration  to 
obtain  a  single  object  estimate  in  which  to  start  the  next  iteration.  This  algorithm 
combines  the  relative  simphcity  of  the  Ayers-Dainty  approach  with  noise  variance 
reduction  through  multiple  image  averaging.  This  modification  to  Ayers-Dainty 
reduces  the  number  of  iterations  required  while  restoring  noisy  speckled  images. 

The  Ayers-Dainty  technique  and  the  modifications  outhned  above  proved  that 
bhnd  deconvolution  was  indeed  possible.  Unfortunately,  they  all  had  a  serious  draw¬ 
back.  Though  they  worked  well  on  simulated  data  where  the  “true”  object  and  PSF 
were  known,  none  of  the  methods  outlined  above  converge  on  the  object  and  PSF  ei¬ 
ther  in  a  mean  squared  error  sense  or  visually  at  consecutive  iterations  (13).  In  other 
words,  if  represents  the  mean  squared  error  at  the  nth  iteration,  then  the  error  at 
a  subsequent  iteration  is  not  necessarily  less  than  the  current  error.  Therefore,  the 
following  is  possible, 

>  4,  (2.10) 

where  contains  the  error  at  the  next  iteration.  Two  methods  have  been  devel¬ 
oped  which  converge  on  either  a  maximum  or  minimum  point;  theoretically  providing 
better  estimates  with  increasing  iterations  and  improved  stopping  criteria  for  the  it¬ 
erative  bhnd  deconvolution  process. 

2.5.2  Maximum-likelihood  estimation.  One  method  has  been  suggested 
which  performs  iterative  bhnd  deconvolution  with  improved  convergence  properties. 
The  goal  of  the  maximum-hkehhood  estimation  (MLE)  technique  is  to  find  the  object 
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intensity  and  the  point  spread  function(s)  that  most  Hkely  created  the  measured  data 
(21).  T.  J.  Holmes  recognized  the  need  to  deconvolve  a  signal  without  knowing  the 
impulse  response  for  confocal  fluorescence  microscopic  imagery  (11).  Holmes  sought 
to  use  MLE  as  a  quantitative  optimization  criterion  to  solve  the  bhnd  deconvolution 
problem  of  quantum-limited  incoherent  imagery. 

In  order  to  utihze  the  iterative  method  of  MLE,  Holmes  derives  a  log-hkeHhood 
function  to  include  the  total  photon  count  and  the  current  object  and  PSF  estimate. 
Each  successive  iteration  produces  a  new  estimate  of  both  the  object  and  PSF  which 
monotonicaUy  increases  the  log-Ukehhood  function.  Holmes  proved  that  since  the 
MLE  method  treated  the  object  and  PSF  arrays  as  probabiHty  density  functions, 
each  estimate  was  impHcitly  constrained  to  be  non-negative  and  retain  unit-volume 
(11).  Unfortunately,  the  simple  unconstrained  MLE  process  did  not  produce  any  re¬ 
sults  resembhng  the  true  object  and  PSF.  Holmes  had  to  apply  an  external  constraint 
to  the  iterative  process. 

As  the  first  constraint.  Holmes  required  that  the  impulse  response  estimate  re¬ 
main  radially  symmetric  at  every  iteration  since  the  PSF  for  many  optical  systems  is 
symmetric.  Outstanding  results  were  obtained  using  the  radial  symmetry  constraint 
on  simulated  images  blurred  with  a  circular  PSF  and  corrupted  by  photon  noise.  As 
Holmes  states,  the  drawback  to  the  radial  symmetry  constraint  is  that  the  trivial 
solution  (where  the  PSF  becomes  a  delta  function  and  the  object  estimate  becomes 
the  measured  image)  wiU  eventually  be  reached  given  enough  iterations  (11). 

To  prevent  the  trivial  solution,  Holmes  tried  restricting  the  bandwidth  of  the 
PSF  to  the  cutoff  frequency  of  the  optical  system.  At  each  iteration  of  the  MLE 
method,  this  constraint  extinguishes  all  information  in  the  OTF  estimate  above  a 
cutoff  frequency  determined  by  the  optics,  then  inverse  Fourier  transforms  the  OTF 
to  obtain  the  PSF,  and  sets  aU  negative  values  in  the  PSF  to  zero  before  returning  the 
new  PSF  back  for  the  next  iteration.  This  constraint  produced  outstanding  results 
in  conjunction  with  MLE  for  restoring  photon-hmited  images.  However,  even  though 
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the  log-likelihood  function  increased  monotonically  using  this  method,  a  point  was 
reached  when  the  object  and  PSF  estimates  began  to  diverge  from  the  true  object 
and  PSF  (11).  Thus,  Holmes  method  has  no  guaranteed  stopping  criteria  which 
present  the  best  estimate  of  both  object  and  PSF. 

T.  J.  Schulz  applied  Holmes’  maximum-hkehhood  estimation  method  to  re¬ 
construct  a  series  of  short-exposure  images  taken  through  atmospheric  turbulence 
without  the  aide  of  a  nearby  point  source  or  any  knowledge  of  the  atmospheric  point 
spread  functions  (21).  Schulz  modified  Holmes’  MLE  method  for  photon-hmited  im¬ 
ages  by  extending  the  iterative  process  to  handle  a  series  of  images  of  the  same  object. 
Noting  that  atmospheric  PSFs  have  neither  radial  symmetry  nor  a  unique  band-limit, 
Schulz  suggested  two  different  modifications  to  the  MLE  method  to  overcome  the 
problem  of  convergence  on  the  trivial  solution.  One  method  utihzes  a  penalized 
maximum-hkehhood  function,  while  the  other  method  is  based  on  the  parameter¬ 
ization  of  the  PSFs  by  phase  errors  distributed  over  an  aperture  (21).  ExceUent 
results  were  obtained  using  simulated  and  actual  photon-hmited  data.  However,  no 
information  was  given  concerning  the  stopping  criteria  for  the  algorithm  or  whether 
the  object  and  PSF  estimates  reached  an  optimal  point  then  diverged  from  the  true 
object  and  PSF. 

2.5.3  Unconstrained  Minimization.  R.  G.  Lane  addressed  the  bhnd  decon¬ 
volution  problem  another  way.  Lane  noted  that  the  Ayers-Dainty  approach  “lacked 
stable  convergence  properties.”  By  minimizing  an  error  metric  which  included  both 
the  object  and  PSF  estimate,  he  modeled  bhnd  deconvolution  as  an  unconstrained 
minimization  problem.  Based  on  a  steepest-descent  search,  Lane’s  algorithm  handles 
additive  noise  and  has  weU-defined  stopping  criteria  (13). 

Lane’s  method  defines  an  error  metric  which  quantifies  how  much  the  current 
estimate  violates  known  constraints.  The  constraints  apphed  in  this  method  are  sim¬ 
ilar  to  those  used  in  previously  mentioned  bhnd  deconvolution  techniques — namely 
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positivity  and  convolution.  Additionally,  Lane  utilizes  a  support  constraint  in  image 
space  which  restricts  the  extent  of  the  object  and  PSF  estimates.  The  combined 
error  function  minimized  in  this  method  is  shown  below, 

Ec  =  Ei  +  Ef,  (2-11) 

where  Ei  and  Ef  represent  the  errors  in  the  image  domain  and  Fourier  domain,  re¬ 
spectively.  With  7/  and  7;,  denoting  the  points  where  f{x,y)  and  h(x,  y)  respectively 
violate  their  image  plane  constraints.  Lane  defines  Ei  as 

=  ff  +  [[  \h{x,y)\^dxdy.  (2.12) 

Ef  quantifies  the  error  in  the  convolution  of  the  two  estimates  with  respect  to  the 
measured  image.  Ef  given  by 

Ef  =  J  J  \G(u,v)  —  F{u,v)H(u,v)\^  dudv,  (2.13) 

where  F{u,v)  and  H{u,v)  are  the  Fourier  transforms  of  the  object  and  PSF  esti¬ 
mates,  respectively  (13). 

Minimization  of  the  combined  error  function  is  accompHshed  through  a  conju¬ 
gate  gradient  minimization  routine.  The  routine  requires  no  additional  constraints 
apart  from  the  objective  function  itself.  Thus,  the  Lane  method  models  an  uncon¬ 
strained  minimization  problem  where  the  pixels  in  the  object  and  PSF  represent 
the  independent  variables  of  the  objective  function  to  be  minimized.  The  fact  that 
this  technique  does  not  use  an  inverse  filter  to  calculate  either  the  object  or  PSF 
estimate  is  one  of  its  strong  points,  since  it  is  not  overly  sensitive  when  either  of  the 
convolutional  components  is  close  to  zero.  Another  advantage  of  conjugate  gradient 
minimization  is  that  the  objective  function  value  monotonically  decreases  with  each 
iteration  until  a  local  minimum  is  found. 
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The  addition  of  a  support  constraint  serves  to  prevent  the  occurrence  of  the 
trivial  solution  where  f{x,y)  is  exactly  g{x,y)  and  the  PSF,  h{x,y),  is  an  impulse 
function,  ^(0, 0).  Constraining  the  object  estimate  to  have  non-zero  values  in  a  region 
smaller  than  the  array  size  ensures  that  the  trivial  solution  has  some  error  associated 
with  it.  Noise  in  the  imaging  process  can  be  modeled  in  the  following  manner, 

9{x,y)  =  f{x,y)  *  h{x,y)  +  c(a;,y),  (2.14) 

where  the  term,  c(a:,y),  represents  additive  noise.  Without  the  support  constraint, 
the  trivial  solution, 

5(a:,2/)  =  sr(a:,y)*5(0,0),  (2.15) 

represents  the  global  minimum  to  Equation  2.14  above. 

Lane  obtained  excellent  results  using  his  method  in  the  presence  of  Gaussian 
and  Poisson  noise  for  simulated  images  using  the  unconstrained  minimization  method 
from  several  different  initial  estimates.  One  problem  the  Lane  method  has  is  the 
relatively  high  probability  of  finding  a  local  minimum  instead  of  the  global  minimum 
depending  on  the  initial  estimates  (13). 

S.  M.  Jefferies  and  J.  C.  Christou  modified  the  Lane  algorithm  (12).  The 
Jefferies  and  Christou  model  includes  not  only  the  positivity,  convolution  and  sup¬ 
port  constraints  utilized  by  previous  blind  deconvolution  algorithms  but  also  applies 
band-limit,  multiple  image  and  Fourier  modulus  constraints.  Jefferies  and  Christou 
claim  their  unconstrained  minimization  technique  achieves  “super-resolution”  for 
continuous  gray-scale  images  and  low  signal-to-noise  ratio  images  by  incorporating 
as  much  a  priori  information  as  possible  in  their  algorithm. 

The  modification  of  the  Lane  algorithm  seeks  to  minimize  an  error  function 
with  four  terms  instead  of  two.  The  error  function,  e,  is  defined  as 

c  =  Ei-\-  Ef  -\-  Ebi  +  Epm,  (2-16) 
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where  the  four  terms  each  represent  an  error  defined  by  the  current  estimate  of  the 
object  and  PSF.  The  image  domain  error  is  exactly  the  same  as  Equation  2.12  in 
the  Lane  method.  The  convolution  error, 

E/ =  —  J J  \G{u,v)  —  F{u,v)H{u,v)f  B(u,v)dudv,  (2-17) 

appears  similar  to  Equation  2.13  with  one  additional  term  containing  a  low  pass  filter, 
B{u,v),  which  ehminates  the  high-spatial  frequency  content  of  the  noise  present  in 
the  measurement.  Thus,  the  object  and  PSF  are  loosely  constrained  to  convolving 
to  a  low  pass  filtered  version  of  the  measured  image,  since  no  information  is  passed 
by  the  pupil  for  spatial  frequencies  greater  than  the  cutoff  frequency,  D/\,  where  D 
is  the  telescope  diameter  and  A  is  the  optical  wavelength. 

The  band-Hmit  error  term,  Ebi,  is  utilized  for  the  first  few  iterations  when 
the  trivial  solution  continually  appears  using  only  the  first  two  error  terms,  then 
convergence  is  obtained  using  the  image  domain  and  convolution  errors.  The  band- 
Hmit  error, 

contains  a  high  pass  filter,  B'{u,v),  which  acts  to  sum  all  frequency  information  in 
the  OTF  outside  the  cutoff  frequency.  Minimizing  the  high-frequency  information 
in  the  PSF  estimate  ensures  that  the  trivial  solution  wiU  not  be  attained.  Since  this 
error  metric  forces  high  spatial  frequency  noise  into  the  object  estimate,  it  has  to  be 
relaxed  before  convergence. 

The  Fourier  modulus  error,  Epm,  utiHzes  Fourier  modulus  information  obtained 
through  speckle  interferometry  to  constrain  the  estimate  of  the  object  spectrum  to 
very  near  the  true  object  spectrum  modulus.  As  expected,  this  constraint  proves 
to  be  a  powerful  tool  in  extracting  an  estimate  very  close  to  the  true  object  from 
noisy  data  as  long  as  the  Fourier  modulus  information  is  available.  When  the  Fourier 
modulus  error  is  used,  the  noise  in  the  convolution  images  mainly  manifests  itself  in 
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the  reconstructed  PSF  leaving  an  object  estimate  almost  void  of  noise.  See  Jefferies 
and  Christou  (12)  for  more  detailed  information  regarding  the  Fourier  modulus  error. 

Few  results  were  obtained  using  the  modified  Lane  method  without  the  Fourier 
modulus  error.  Those  shown  were  compared  against  the  excellent  estimates  made 
with  the  Fourier  modulus  error.  Jefferies  and  Christou  performed  their  analysis  on 
images  corrupted  by  zero-mean  Gaussian  noise.  What  remains  to  be  asked  are  the 
following: 

1.  Can  the  addition  of  band-limit  information  to  Lane’s  blind  deconvolution 
technique  perform  well  in  the  presence  of  photon  noise? 

2.  Can  this  method  be  effectively  used  on  compensated  imagery  from  an  adap¬ 
tive  optics  system? 

3.  Since  neither  Lane  nor  Jefferies  and  Christou  discuss  the  support  constraint 
in  detail,  how  does  varying  the  support  region  affect  the  restoration  of  images? 

2. 6  Summary 

This  chapter  reviewed  fundamentals  of  imaging  through  atmospheric  turbu¬ 
lence  and  presented  a  brief  explanation  of  adaptive  optics.  Since  imagery  from 
adaptive  optics  stiU  requires  post-processing,  the  method  of  blind  deconvolution  was 
reviewed  explicitly.  Blind  deconvolution  involves  restoration  of  an  image  or  a  series  of 
images  without  prior  knowledge  of  the  blurring  function  or  PSF.  The  Ayers-Dainty 
approach  was  the  first  workable  algorithm  toward  solving  the  bUnd  deconvolution 
problem.  Due  to  the  inabihty  of  Ayers-Dainty  to  converge  on  a  “best”  estimate, 
two  other  methods  have  been  proposed.  Maximum-likeUhood  estimation  uses  an 
iterative  method  which  increases  the  likelihood  that  the  object  and  PSF  correctly 
estimate  the  “true”  object  and  PSF.  Results  from  Holmes,  however,  reveal  that  the 
maximum-likelihood  estimation  reaches  an  optimal  estimate  then  diverges  from  the 
“true”  object  and  PSF  with  further  iterations.  Lane  proposed  an  unconstrained 
minimization  solution  to  the  blind  deconvolution  problem.  Constraints  are  effec- 
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tively  modeled  in  an  error  function  which  can  be  minimized  through  the  conjugate 
gradient  method.  Further  modifications  to  the  Lane  method  developed  by  Jefferies 
and  Christou  are  worthy  of  continued  analysis.  Chapter  3  discusses  the  modified 
Lane  technique  for  experimentation  on  the  effects  of  support  size  and  photon  noise 
on  image  restoration. 
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III.  Methodology 


3. 1  Introduction 

This  chapter  transitions  from  a  basic  understanding  of  the  unconstrained  min¬ 
imization  technique  for  bhnd  deconvolution  to  its  implementation  on  simulated  im¬ 
ages.  In  order  to  fuUy  understand  the  approach  taken,  the  notation  for  several  terms 
are  defined  followed  by  a  defense  for  computer  simulation  of  this  problem.  The  latter 
part  of  the  chapter  includes  several  key  factors  in  the  development  of  experimenta¬ 
tion  to  analyze  the  Jefferies  and  Christou  modification  to  the  Lane  method  called 
the  iterative  deconvolution  algorithm  (IDA).  First,  the  goals  of  subsequent  experi¬ 
mentation  and  testing  are  listed.  Then,  a  discussion  of  how  the  functions  affecting 
the  bhnd  deconvolution  algorithm  is  presented  followed  by  a  discussion  of  an  error 
metric  utihzed  in  the  comparison  of  different  results.  Finally,  a  fisting  of  the  vari¬ 
ables  affecting  the  blind  deconvolution  process  emphasizes  which  ones  will  be  altered 
and  which  ones  will  remain  constant  for  the  purposes  of  experimentation. 

3.2  Terminology 

As  alluded  to  in  the  previous  chapter,  imaging  is  modeled  as  a  linear,  shift 
invariant  process  with  the  following  form: 

9{x,y)  =  f{x,y)*h{x,y),  (3.1) 

where  {x,y)  represents  an  arbitrary  point  in  the  image  domain,  denotes  the  two- 
dimensional  convolution  of  the  “true”  object,  f{x,y),  and  the  point  spread  function 
(PSF),  h{x,y),  and  g{x,y)  defines  the  measured  image.  Departures  from  the  ideal 
convolution  include  noise  inherent  in  the  imaging  process  and  assorted  non-linearities 
(2).  For  short-exposure  images,  photon  noise  manifests  itself  in  the  random  arrival 
of  photons  at  the  detector.  AU  of  these  departures  can  be  added  to  Equation  3.1  in 
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the  form  of  a  single  contamination  term.  A  new  mathematical  model  of  the  imaging 
process  follows, 

9{x,y)  =  fix,y)  *  h{x,y)  +  c{x,y),  (3.2) 

where  c{x,  y)  represents  the  contamination  term.  The  problem  of  blind  deconvolu¬ 
tion  is  to  determine  f{x,y)  given  only  the  measured  image,  g{x,y),  and  very  Httle 
information  about  h[x,y)  and  c(x,y).  Subsequent  sections  discuss  each  of  the  terms 
in  the  imaging  process  and  how  they  relate  to  IDA  using  the  unconstrained  mini¬ 
mization  approach. 

3.2.1  Object.  Obtaining  a  close  representation  of  the  “object”  is  the  goal 
of  any  imaging  process,  though  the  object  itself  can  never  be  fully  reahzed.  Even 
when  one  is  close  enough  to  observe  an  object  with  the  unaided  eye,  the  optical 
quaUties  of  the  human  visual  system,  like  the  finite  extent  of  the  pupil,  distort  the 
actual  object.  For  the  purpose  of  analysis  here,  the  “true”  object,  f{x,y),  is  a  given 
array  of  irradiance  values.  Simulated  blurred  images  are  created  using  Equation  3.2 
above,  and  the  results  from  IDA  are  analyzed  for  comparison  purposes  using  the 
true  object.  The  term  “object  estimate,”  hereafter  designated  f{x,y),  represents 
the  best  guess  at  the  “true”  object  that  IDA  can  produce.  The  Fourier  transforms  of 
both  the  true  object  and  the  object  estimate  are  represented  by  the  object  spectra, 
F[u,v)  and  F{u,v),  respectively,  where  {u,v)  represents  an  arbitrary  point  in  the 
spatial  frequency  domain. 

3.2.2  Point  Spread  Function  (PSF).  Sometimes  referred  to  as  the  impulse 
response,  the  true  PSF,  represented  by  h{x,y),  contains  the  blurring  function  for 
the  entire  optical  system.  This  includes  blurring  caused  by  atmospheric  turbulence 
and  the  optics.  For  adaptive  optics  compensated  images,  the  PSF  includes  blurring 
from  the  optics  and  any  residual  phase  errors  uncorrected  by  the  adaptive  optics 
system  (21).  IDA  estimates  the  PSF  which  wiU  be  denoted  as  h[x,y).  The  Fourier 
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transform  of  the  PSF  is  called  the  optical  transfer  function  (OTF).  The  true  and 
estimated  OTFs  are  represented  by  H{u,v)  and  H{u,v),  respectively. 

3.2.3  Image.  The  image  is  equivalently  referred  to  as  the  convolution 
image,  the  measured  image  or  the  measurement.  The  function,  g[x,y),  represents 
the  magnitude  of  the  convolution  plus  any  contamination  as  revealed  by  Equation 
3.2.  Its  Fourier  transform  is  called  the  image  spectrum  and  is  denoted  by  G{u,v). 

3.2.4  Positivity.  Positivity  is  the  assumption  that  both  the  object,  f{x,y), 
and  the  PSF,  h(x,y),  have  positive  intensity  values.  One  term  in  the  algorithm’s 
objective  function  constrains  the  object  and  PSF  estimates  to  have  a  minimum  of 
negative  intensity  values. 

3.2.5  Support.  Support  is  the  only  “hard”  constraint  utilized  by  the  IDA 
method.  Two  arrays,  representing  the  extent  of  the  object  and  PSF,  respectively, 
must  be  provided  for  the  algorithm  along  with  the  measured  image  and  initial  esti¬ 
mates  of  the  object  and  PSF.  These  arrays  constrain  the  object  and  PSF  estimates  to 
have  non-zero  values  only  where  the  support  array  has  a  value  of  one.  The  remainder 
of  the  image  plane  is  set  to  zero  corresponding  to  the  complementary  region  where 
the  support  array  is  zero.  Since  the  pixels  within  the  object  and  PSF  estimates  be¬ 
come  the  variables  for  minimization  in  the  unconstrained  minimization  process,  the 
support  arrays  actually  reduce  the  number  of  variables  IDA  is  required  to  minimize. 

3.2.6  Error  Metric.  The  error  metric  is  a  function  which  measures  how 
close  or  how  far  an  estimate  is  in  relation  to  its  true  value.  A  separate  error  met¬ 
ric  will  be  calculated  for  each  object  and  each  PSF  estimate.  Since  this  metric  is 
calculated  using  the  true  object  or  PSF,  it  will  not  be  useful  for  processing  real  im¬ 
agery.  However,  using  Davey’s  philosophy,  such  an  error  metric  serves  as  an  excellent 
preliminary  assessment  of  the  algorithm’s  performance  (4). 
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3.3  Computer  Simulation  of  the  Problem 

This  investigation  utilizes  simulation  to  reconstruct  computer  generated  images 
for  several  reasons.  First,  baseline  data  is  necessary  using  pristine,  infinite  signal 
images  to  determine  a  maximum  ability  to  reconstruct  an  object  under  a  given 
set  of  conditions.  Second,  the  effects  of  photon  noise  are  to  be  studied  in  detail. 
Obtaining  actual  imagery  where  photon  noise  is  dominant  is  possible,  but  only  with 
simulated  images  can  the  photon  level  be  carefully  varied.  Modeling  photon  arrivals 
as  a  Poisson  process  in  a  computer  simulation  allows  one  to  analyze  the  performance 
of  IDA  under  the  strict  condition  where  the  contamination,  c[x,y)  in  Equation  3.2, 
contains  only  contributions  due  to  photon  noise.  Finally,  computer  generated  data 
can  be  processed  and  analyzed  on  the  same  machine  which  means  that  the  complexity 
of  data  transfer  is  avoided. 

Experimentation  through  computer  simulation  is  similar  to  laboratory  experi¬ 
mentation.  When  one  approaches  computer  simulation  using  the  scientific  method, 
tangible  results  can  be  obtained.  Therefore,  before  jumping  head-first  into  the  itera¬ 
tive  deconvolution  algorithm,  several  questions  must  be  asked  and  answered  to  guide 
the  experimental  process.  Five  questions  follow  which  provide  a  scientific  approach 
to  computer  simulation  testing  (23). 

1.  What  is  the  experimental  objective? 

2.  What  phenomena  will  be  observed? 

3.  How  will  observations  be  quantified? 

4.  What  are  the  independent  variables  affecting  the  phenomena? 

5.  What  factors  wiU  be  held  constant  in  the  experiments? 

Subsequent  discussion  seeks  to  answer  each  of  the  above  questions.  The  main 
purpose  of  the  next  several  sections  is  to  give  the  goals  for  the  experimentation 
discussed  in  Chapter  4.  In  addition,  key  factors  and  variables  in  the  unconstrained 
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minimization  process  are  listed  and  defined  to  provide  a  better  understanding  of  the 
results  obtained  as  well  as  the  assumptions  that  apply  to  those  results. 

3.4  Experimental  Objectives 

The  overall  driving  concern  for  this  study  lies  in  answering  the  following  ques¬ 
tion:  Can  unconstrained  minimization  eflFectively  deconvolve  an  unknown  object  and 
unknown  point  spread  function  from  a  single  measurement?  Though  a  simple  “yes” 
or  “no”  could  answer  the  question,  neither  question  nor  answer  is  simple  in  this 
case.  Implementation  of  unconstrained  minimization  through  the  IDA  technique  al¬ 
lows  many  different  aspects  of  the  reconstruction  process  to  be  estimated  or  varied. 
Therefore,  several  variables  were  chosen  for  further  analysis.  The  objectives  relating 
these  variables  are  hsted  below. 

1.  Develop  a  baseline  for  object  estimates  using  infinite  signal  measurements, 
and  determine  the  optimal  support  size  necessary  for  reconstruction.  This  first  ob¬ 
jective  will  allow  IDA  to  process  several  noiseless,  blurred  images  using  a  series  of 
well-defined  support  constraints  so  that  generalizations  can  be  made  across  all  the 
images  tested. 

2.  Determine  the  effect  of  photon  noise  on  the  object  estimate  and  compare 
with  results  previously  obtained  by  Lane.  Photon  noise  can  be  considered  the  origin 
of  shot  noise  in  detectors  that  emit  an  electron  upon  absorption  each  photon  (3). 
Lane  revealed  that  decreasing  the  light  level  present  in  the  measurement  significantly 
affects  the  reconstruction  (13).  This  second  objective  seeks  to  test  the  modifications 
made  by  Jefferies  and  Christou,  specifically  the  addition  of  a  band-limiting  mask 
in  the  convolution  error  shown  by  Equation  2.16  (12).  Results  using  this  band- 
limit  mask  can  be  compared  with  results  obtained  using  Lane’s  method  without  the 
band-limit  mask. 

3.  Determine  the  utility  of  the  band-limit  error  term  toward  improving  the 
reconstruction  process.  Jefferies  and  Christou  define  the  third  error  term  in  their  ob- 
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jective  function  as  the  band-limit  error  which  minimizes  high  frequencies  in  the  OTF 
estimates  (12).  This  third  objective  seeks  to  analyze  whether  use  of  the  band-hmit 
error  shown  by  Equation  2.18  contributes  in  a  positive  way  toward  the  reconstruction 
process  either  in  terms  of  reduced  iterations  or  improved  object  and  PSF  estimates. 

The  enumerated  objectives  just  listed  form  the  main  goals  of  the  three  ex¬ 
periments  for  which  results  are  presented  and  analyzed  in  the  next  chapter.  The 
purpose  for  those  experiments  will  be  to  analyze  the  performance  of  Jefferies  and 
Christou’s  IDA  technique.  To  better  understand  the  results,  however,  one  first  must 
comprehend  the  algorithm  itself  and  how  it  achieves  a  best  estimate  of  an  object  and 
PSF.  Thus,  time  must  be  devoted  to  understanding  the  unconstrained  minimization 
technique  for  blind  deconvolution. 

8.5  Constraint  Driven  Phenomena 

Laboratory  experiments  are  usually  designed  to  test  or  reveal  some  physical 
phenomenon.  Computer  simulation  can  accomphsh  a  similar  goal  if  the  mathemat¬ 
ical  model  represents  a  close  approximation  to  reality.  Therefore,  the  phenomena 
revealed  in  the  subsequent  experiments  result  directly  from  the  computer  algorithm. 
Although  blind  deconvolution  was  defined  and  a  summary  of  different  techniques  was 
outhned  in  Chapter  2,  a  more  complete  knowledge  of  the  unconstrained  minimization 
technique  utihzed  in  this  investigation  is  necessary.  IDA,  developed  by  Jefferies  and 
Christou  as  a  modification  to  Lane’s  method,  implements  this  proposed  bhnd  decon¬ 
volution  solution  and  is  the  object  for  subsequent  analysis.  The  following  paragraphs 
seek  to  give  insight  into  the  Jefferies  and  Christou  algorithm. 

Since  the  Ayers-Dainty  method  had  no  well-defined  stopping  criteria,  Lane 
sought  an  unconstrained  minimization  approach  to  bhnd  deconvolution  by  devel¬ 
oping  an  algorithm  which  converged  on  a  best  set  of  object  and  PSF  estimates. 
The  Lane  technique  formulates  the  same  constraints  appHed  at  each  iteration  in  the 
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Ayers-Dainty  method  into  a  single  error  function.  This  function  is  then  minimized 
using  conjugate  gradient  minimization. 

3.5.1  Minimization.  Conjugate  gradient  minimization  (CGM)  is  one  of 
a  number  of  steepest-descent  iterative  search  methods.  CGM  finds  at  least  a  local 
minimum  for  a  multi- variate  objective  function  given  an  initial  starting  point  for  each 
variable.  With  a  first  guess,  the  CGM  method  calculates  the  value  of  the  objective 
function  as  well  as  the  gradient  at  that  location.  Using  this  information,  the  method 
then  selects  a  subsequent  guess  in  a  direction  opposite  to  the  gradient.  Thus,  the  next 
guess  traverses  the  steepest-descent  along  the  complex,  multidimensional  function. 
With  this  new  guess,  a  new  objective  function  value  and  gradient  are  found.  This 
process  is  repeated  until  the  new  objective  function  value  is  greater  than  the  previous, 
which  implies  that  the  previous  point  was  a  local  minimum. 

IDA  implements  conjugate  gradient  minimization  on  a  single  error  function 
with  several  terms.  The  pixel  values  within  the  support  of  the  object  and  PSF 
estimate  arrays  become  the  variables  to  be  minimized.  Therefore,  CGM  is  performed 
on  a  single  objective  function  containing  hundreds  to  thousands  of  variables  for  even 
the  relatively  small  image  plane  array  of  64  x  64.  Fortunately,  the  error  terms  derived 
by  Lane  and  Jefferies  have  analytic  derivatives  which  simplify  the  calculation  of  the 
gradient  at  each  iteration.  Though  understanding  the  CGM  process  is  important, 
the  error  terms  are  the  key  to  the  unconstrained  minimization  process. 

3.5.2  Image  Error.  As  discussed  in  Section  2.5.3,  the  image  domain  error 

term, 

=  / /  \fix,y)\^dxdy  +  f f  \h{x,y)\^dxdy,  (3.3) 

quantifies  the  squared  magnitude  of  the  negative  pixels  inside  the  support  regions  for 
both  the  object  and  PSF  where  7/  and  jh  respectively  represent  the  set  of  negative 
pixels  within  the  object  and  PSF  support  regions.  Instead  of  forcing  all  negative 


3-7 


values  to  zero  as  is  done  in  the  Ayers-Dainty  approach,  IDA  seeks  to  find  estimates  of 
the  object  and  PSF  where  the  magnitudes  of  the  negative  pixel  values  are  minimized. 

3.5.3  Convolution  Error.  The  convolution  error  term, 

F?/ = -^  J J  \G(u,v)  —  F(u,v)H{u,v)\^  B{u,v)  dudv,  (3-4) 

where  represents  the  total  number  of  pixels  in  the  image  plane  and  B{u,v)  is 
a  band-Hmit  mask,  should  always  dominate  the  combined  error  of  the  objective 
function  since  the  deconvolution  of  the  object  and  PSF  from  the  measured  image  is 
the  goal  of  this  process.  This  error  term  can  best  be  expressed  in  terms  of  the  Fourier 
transforms  due  to  the  convolution  theorem  (7).  As  such.  Equation  3.4  is  minimized 
when  the  product  of  both  F{u,v)  and  H{u,v)  equals  G{u,v).  With  the  addition  of 
noise  in  the  imaging  process,  G{u,v)  contains  high  frequency  values  not  present  in 
either  F{u,v)  or  H{u,v).  Thus,  a  band-limit  mask,  B{u,v),  is  added  to  prevent  the 
product  of  F{u,v)  and  H{u,v)  from  having  to  reconstruct  the  high  frequency  noise 
content.  This  band-Hmit  mask  should  allow  IDA  to  handle  a  greater  amount  of  noise 
than  the  Lane  method,  which  does  not  include  the  mask. 

Another  point  to  be  considered  is  that  the  noise  present  within  the  passband 
will  have  to  be  accounted  for  within  the  object  and  PSF  estimates.  Even  the  best 
estimate  of  the  object  and  PSF  will  contain  a  certain  amount  of  error  resulting 
from  the  noise  whose  spectral  content  Hes  within  the  cutoff  frequency.  Thus,  the 
fundamental  Hmit  for  restoring  the  measured  image  is  not  reached  until  the  signal 
cannot  be  retrieved  from  the  noise  present  within  the  passband. 

3.5.4  Band-Limit  Error.  The  band-Hmit  error  term, 

^61  =  -^  //  (3-5) 


3-8 


is  minimized  when  the  magnitude  of  all  spatial  frequency  values  in  the  OTF  beyond 
1.39  times  the  cutoff  frequency  is  zero  or  very  small,  where  B'{u,v)  represents  a 
high  pass  filter  which  blocks  out  OTF  values  within  1.39  times  the  cutoff  frequency 
(12).  This  error  term  eliminates  the  occurrence  of  the  trivial  solution  by  preventing 
the  PSF  estimate  from  approaching  a  delta  function.  Since  the  delta  function  has  a 
Fourier  transform  extending  out  to  infinite  spatial  frequencies,  hmiting  the  OTF  to 
a  specific  passband  rejects  the  trivial  solution. 

Use  of  this  error  term  also  has  a  distinct  disadvantage  which  was  alluded  to  by 
Jefferies  and  Christou  (12).  As  stated  above  in  Section  3.5.3,  the  spectral  content  of 
the  noise  within  the  passband  must  manifest  itself  in  the  object  and  PSF  estimates. 
When  the  PSF  estimate  is  constrained  by  the  band-hmit  error  term,  almost  all  of  the 
noise  within  the  passband  is  accounted  for  in  the  object  estimate.  This  is  the  reason 
why  Jefferies  states  that  this  constraint  “has  to  be  relaxed  before  convergence”  (12). 

3.6  Defining  Comparison  Metrics 

Metrics  are  necessary  to  quantify  errors  in  the  reconstructed  objects  and  PSFs 
and  for  comparison  against  other  object  and  PSF  estimates.  Since  this  investigation 
solely  uses  computer  generated  images  where  the  “true”  object  and  PSF  are  known 
exactly,  metrics  are  presented  which  compare  an  estimate  of  either  the  object  or  PSF 
with  its  actual  value. 

Goodman  and  Belsher  [1976]  discuss  an  error  metric  which  uses  a  filtered  ob¬ 
ject,  call  it  f{x,y),  rather  than  the  true  object,  f{x,y),  since  “restoration  of  an 
object’s  frequency  components  beyond  the  diffraction  limited  cutoff  of  the  optical 
system  is  impossible  to  achieve  with  any  linear  invariant  restoration  filter”  (9).  Thus, 
Goodman  suggests  that  a  low  pass  filtered  version  of  the  true  object  is  the  best  that 
can  be  obtained  during  post-processing.  With  the  understanding  that  only  those 
frequency  components  within  the  passband  can  be  restored,  a  new  error  metric  for 
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the  object,  follows: 


_  /Xy  \K{U,V)  -  Fn{u,v)\^dudv 

SJ^\Fr.{u,v)\^dudv  ’ 

where  7  includes  all  points  within  the  diffraction  hmited  cutoff,  Fn{u,v)  represents 
the  normahzed  spectrum  of  the  filtered  true  object  and  Fn{u,v)  contains  the  nor- 
mahzed  spectrum  of  the  object  estimate.  Similarly,  the  PSF  error  can  be  written  as 
a  function  of  the  true  OTF,  H(u,v),  and  the  normalized  OTF  estimate,  Hn{u,v), 
shown  by 

>  ^  Il-y  \H{U,V)  -  Hr,{u,v)\^  dudv 

^  JJ^\H(u,v)\^  dudv 

These  error  metrics  count  as  error  only  those  differences  between  the  restored  spec¬ 
trum  and  the  filtered  true  spectrum  that  lie  within  the  passband.  Division  by  the 
integral  over  the  power  spectrum  expresses  the  comparison  metric  as  a  fraction, 
for  estimates  that  are  reasonably  close  to  the  actual  value. 

Since  aU  objects  used  in  this  investigation  will  contain  point  sources  of  varying 
intensity,  additional  metrics  are  useful  in  defining  the  error  in  point  source  location 
and  the  intensity  ratio  for  binary  stars  in  the  object  estimates.  For  both  true  objects 
and  images,  the  intensity  ratio  for  a  binary  star  is  often  denoted  as  1  :  77,  where  rj  is 
defined  as  the  ratio  of  the  intensity  of  the  dimmer  star  to  the  intensity  of  the  brighter 
star.  Thus,  77,  by  definition,  is  less  than  1. 

A  location  error  is  defined  for  this  investigation  as  the  deviation  in  point  source 
location  between  the  true  object  and  its  estimate  divided  by  the  blur  radius,  $.  $ 
represents  the  radial  distance  from  the  center  of  the  array  to  the  first  zero  in  the 
true  PSF,  when  the  PSF  is  an  Airy  function.  If  the  PSF  remains  constant  for  all 
images,  then  the  following  equation  can  be  used  to  calculate  the  location  error  for  a 
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particular  point  source, 


where  {xp,yp)  represents  the  location  of  point  source,  p,  in  the  “true”  object  and 
(xp,  Pp)  represents  the  location  of  point  source,  p,  in  the  object  estimate. 

3. 7  Independent  Variables 

From  the  objectives  selected  in  Section  3.4,  several  aspects  of  the  unconstrained 
minimization  approach  to  bhnd  deconvolution  need  to  be  isolated.  By  varying  a  sin¬ 
gle  component  within  the  isolated  set,  conclusions  can  be  drawn  regarding  the  effect 
of  that  variable  on  the  reconstruction  process.  To  meet  the  enumerated  objectives, 
this  investigation  will  analyze  the  IDA  technique  by  varying  the  object  contents,  the 
support  size  and  the  light  level  present  in  each  image. 

3. 7. 1  Objects  and  Image  Measurements.  Four  different  objects  were  chosen 
to  study  the  abihty  of  IDA  to  reconstruct  blurred  images  of  astronomical  objects. 
The  four  object  arrays,  64  X  64  pixels  in  size,  contain  point  sources  of  differing  inten¬ 
sity.  Object  A,  shown  in  Appendix  B  as  Figure  B.l,  represents  a  single  star  centered 
in  an  empty  field.  Object  B,  Figure  B.3,  is  a  binary  star  object  of  intensity  ratio  1:0.7 
with  a  very  wide  separation.  Object  C,  Figure  B.5,  shows  another  binary  star  object 
with  intensity  ratio  1:0.7  where  the  separation  is  closer.  The  final  object.  Object  D 
(Figure  B.7),  has  the  same  binary  star  object  with  a  very  close  separation.  Instead 
of  using  separate  PSFs,  binary  point  sources  of  varying  separation  are  employed  to 
show  the  resolving  power  of  the  unconstrained  minimization  technique. 

A  single  PSF  was  selected  to  blur  each  of  the  objects  hsted  above.  This  PSF 
has  a  certain  radius,  $,  from  its  center  to  the  first  zero  which  can  measured  in  pixels, 
micro-radians  or  arc-seconds.  By  defining  the  separation  distance  between  the  point 
sources  in  Objects  B,  C,  and  D  as  0,  then  the  ratio,  represents  a  general  means 
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of  quantifying  the  separation  in  the  three  objects.  Using  the  ratio,  the  separation 
between  the  point  sources  for  Objects  B,  C  and  D  are  respectively,  1.81,  0.84  and 
0.71.  When  the  PSF  is  convolved  with  the  true  point  source  objects,  much  of  the 
detail  disappears  as  revealed  in  Figures  B.2,  B.4,  B.6  and  B.8  which  correspond 
directly  to  Objects  A,  B,  C  and  D. 

5.7.2  Support  Size.  The  first  experiment  tests  the  abihty  of  the  algorithm 
to  reconstruct  noise-free  blurred  images  with  differing  support  sizes.  Support  is  the 
only  firm  constraint  input  into  IDA;  thus  determining  an  optimal  support  region 
in  the  first  experiment  is  necessary  to  effectively  running  other  experiments.  It  is 
desirable  to  test  the  correlation  of  support  region  size  with  definable  areas  within 
the  impulse  response  or  true  PSF.  Since  a  single  PSF  in  the  form  of  an  Airy  function 
is  used  to  blur  each  object,  support  size  will  vary  according  to  the  “zeros”  of  the 
PSF.  Three  different  support  sizes,  given  in  terms  of  the  radius  of  a  circle,  were 
chosen.  The  first  support  size  extends  to  the  zero  region  surrounding  the  main  lobe 
of  the  true  PSF,  having  the  same  radius,  $.  The  second  and  third  support  sizes 
respectively  correspond  to  the  nuU  regions  just  after  the  first  and  second  “rings” 
of  the  PSF.  Support  areas  are  circular  regions  centered  for  the  PSF  support  and 
centered  on  each  point  source  location  for  each  object  support.  Thus  the  object 
support  for  one  of  the  binary  stars  is  the  union  of  two  circles  of  varying  size  each 
centered  on  the  point  source  pixel  locations.  Figure  3.1  reveals  an  example  of  a 
support  region  for  Object  B. 

Experimental  testing  will  utiHze  the  three  separate  support  regions  defined  for 
each  object  as  input  into  the  unconstrained  minimization  method.  Results  from  this 
testing  can  be  analyzed  to  determine  how  support  size  affects  the  bhnd  deconvolution 
reconstruction  process. 

3.7.3  Average  Photo  Events,  K.  The  most  basic  source  of  noise  lies  in 
photon  fluctuations  associated  with  the  detection  of  the  finite  amount  of  fight  energy 
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Figure  3.1  The  object  support  region  above  relates  to  Object  B  where  each  circle 
is  centered  on  one  of  the  two  stars  in  the  true  object. 

available  in  the  imaging  system.  These  photon  fluctuations  pose  a  fundamental 
limitation  to  the  “restorabihty”  of  a  degraded  image  (9).  Photon  arrival  at  the 
detector  has  been  previously  modeled  as  a  Poisson  process  (9,  13,  20).  A  similar 
model  will  be  used  here  to  create  simulated  images  corrupted  by  photon  noise.  To 
understand  the  model,  first  one  must  comprehend  the  factors  affecting  photon  arrival. 
The  mean  number  of  photons  or  the  average  number  of  photo  events  is  a  function 
of  the  visual  magnitude  of  the  object,  the  image  exposure  time,  the  mean  imaging 
wavelength  of  the  detector  and  the  Hght  gathering  capacity  of  the  optical  device  (3). 

The  model  used  to  generate  photon  limited  images  makes  use  of  the  Poisson 
probability  mass  function  defined  by  Ross  [1993], 

p(®)  =  (^)  ,  *  =  0,1,2,...  (3.9) 

where  x  represents  the  number  of  photon  arrivals  at  a  certain  location  and  A  corre¬ 
sponds  to  the  mean  number  of  photons  which  arrive  at  a  detector  cell  over  a  specified 
period  of  time  (20).  The  number  of  photo  events  at  each  pixel  location  is  a  random 
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Table  3.1  Visual  magnitude  and  photo  events  for  common  sky  objects. 


OBJECT 

rriy 

K 

Venus 

-4.4 

2.35 

10" 

Jupiter 

-2.5 

4.11 

10® 

Sirius 

-1.5 

1.59 

10® 

Artificial  SateUite  (typical) 

-.96 

1.00 

10® 

Polaris 

+2.2 

5.54 

lO'* 

variable  from  the  Poisson  distribution  where  the  ratio  of  the  irradiance  in  a  pixel 
to  the  total  energy  in  the  noise-free  image  is  proportional  to  the  mean  of  the  ran¬ 
dom  variable,  A.  The  proportionality  constant,  K,  represents  the  average  number  of 
photo  events  in  the  entire  image  over  a  single  exposure  or  integration  time.  Thus,  the 
mean  of  the  Poisson  random  variable  at  each  pixel  location,  A(a:,t/),  can  be  expressed 
in  the  following  way. 


\{x,y)  = 


Kg{x,y) 


(3.10) 


IS  9{^,y)dxdy' 

where  g{x,y)  represents  a  point  in  the  noise-free  convolution  of  f{x,y)  and  h{x,y). 
Equation  3.10  above  was  used  to  generate  images  corrupted  by  photon  noise  for 
differing  photo  event  levels,  K.  Apparent  visual  magnitude,  to„,  is  a  convention 
employed  by  astronomers  to  compare  the  relative  brightness  of  objects  in  the  night 
sky  (17).  Each  unit  decrease  in  visual  magnitude  corresponds  to  a  2.5  factor  increase 
in  visual  brightness;  thus,  smaller  indicates  a  brighter  object.  This  same  visual 
magnitude  can  be  used  to  calculate  the  average  photon  flux  of  an  object.  Table  3.1 
reveals  the  number  of  photo  events  corresponding  to  the  visual  magnitude  for  several 
common  sky  objects  (23),  based  on  a  1  meter  aperture,  an  imaging  wavelength  of 
500  nm  and  a  1.8-ms  integration  time.  Three  different  values  were  selected  for  K: 
100,000,  10,000,  and  5,000.  The  photon  hmited  images  for  the  four  objects  appear 
in  Appendix  B,  Figures  B.IO  through  B.21. 
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3. 8  Constants 


Now  that  several  variables  have  been  reviewed,  a  couple  of  factors  must  be 
held  constant  for  accurate  analysis  of  the  independent  variables.  The  constants  for 
subsequent  experimentation  are  the  OTF  and  the  initial  estimates  input  into  Jefferies 
and  Christou’s  bhnd  deconvolution  algorithm. 

3.8.1  Optical  Transfer  Function.  As  already  stated,  the  OTF  used  to  create 
all  simulated  images  is  held  constant  so  that  comparisons  can  be  made  regarding  the 
abihty  of  IDA  to  reconstruct  it  for  different  objects  and  differing  amounts  of  noise. 
Subsequent  experiments  use  the  Airy  pattern  for  the  true  PSF  shown  in  Figure 
3.2.  This  PSF  is  theoretically  the  PSF  of  a  circular  diffraction-hmited  system  (11). 
Goodman  [1968]  defines  the  Fourier  transform  of  this  PSF,  or  rather  the  OTF,  for  a 
diffraction-hmited  circular  aperture  for  incoherent  hght  imaging  as  follows: 


2 

cos  ^  (f-)  -  f-\/l  -  (f-) 

a p  <po 

H(p)  =  . 

IT 

\Po  /  po  y  \Po  / 

(3.11) 

0 

otherwise. 

where  p  =  y/u^  +  and  po  denotes  the  diffraction-hmited  cutoff  frequency  (7).  By 
keeping  this  factor  constant,  any  deviations  observed  in  the  PSF  estimate  result 
primarily  from  the  object  or  the  noise. 

3.8.2  Initial  Guesses.  The  iterative  deconvolution  algorithm  requires  an 
initial  estimate  at  both  the  object  and  PSF  for  a  given  measured  image.  In  the  past, 
others  (1,  12,  13)  have  input  white  noise  for  initial  estimates  of  both  the  object  and 
PSF.  AU  runs  in  this  investigation  wiU  use  a  single  Gaussian  estimate  for  the  PSF 
and  wiU  use  the  measurement  as  the  initial  object  estimate.  Though  this  procedure 
draws  close  to  the  trivial  solution,  it  also  represents  a  more  logical  first  guess  at  both 
values.  The  object  should  simply  be  a  variation  on  the  measured  image  whereas  the 
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Figure  3.2  The  PSF  shown  above  was  convolved  with  the  four  true  objects  to  obtain 
the  blurred  images  found  in  Appendix  B. 

PSF  should  contain  some  of  the  same  properties  as  the  Gaussian — radial  symmetry 
and  a  smooth  roU-off. 

3.9  Summary 

The  merits  of  the  unconstrained  minimization  approach  to  bhnd  deconvolution 
derived  by  Lane  and  modified  by  Jefferies  and  Christou  warrants  further  analysis. 
This  chapter  discussed  several  points  crucial  to  understanding  the  experiments  for 
this  investigation.  As  presented,  the  terminology  listed  requires  merely  a  basic  un¬ 
derstanding  of  Fourier  optics.  The  simplifying  assumptions  made  and  the  models 
developed  here  allow  direct  application  of  the  findings  of  this  study  to  actual  astro¬ 
nomical  images.  The  objectives  defined  in  this  chapter  are  further  developed  into 
actual  experiments  in  Chapter  4.  Additionally,  the  following  chapter  presents  the 
results  of  those  experiments  and  analyzes  the  results  in  light  of  the  terminology, 
assumptions  and  metrics  just  discussed. 
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IV.  Results  and  Analysis 


4.1  Introduction 

To  determine  how  well  Jefferies  and  Christou’s  iterative  deconvolution  algo¬ 
rithm  (IDA)  can  reconstruct  compensated  astronomical  imagery,  several  experiments 
are  performed.  These  experiments  utilize  computer  simulation  to  analyze  the  algo¬ 
rithm’s  performance  under  specific  controlled  conditions.  Post-processing  is  per¬ 
formed  using  computer  generated  images.  The  following  sections  provide  detailed 
information  regarding  the  three  experiments  and  draw  conclusions  from  the  results 
of  each  experiment. 

4.2  Experiment  One 

Develop  a  baseline  for  object  estimates  using  infinite  signal  measurements,  and 
determine  the  optimal  support  size  necessary  for  reconstruction.  As  stated  in  Section 
3.4,  the  first  major  objective  of  this  experiment  involves  testing  IDA  on  noise-free 
images  to  understand  the  capability  of  this  bhnd  deconvolution  technique  without 
the  presence  of  noise.  A  secondary  goal  involves  varying  the  support  constraint  to 
determine  its  impact  on  the  reconstruction  process. 

4.2.1  Simulation  Parameters.  Input  images  for  this  experiment  represent 
the  convolution  of  the  “true”  objects  with  a  single  PSF  as  discussed  in  Section 
3.6.1.  These  four  noise-free  images  (Figures  B.2,  B.4,  B.6  and  B.8)  are  displayed  in 
a  three-dimensional  view  where  increasing  irradiance  corresponds  to  a  higher  value 
on  the  graph.  In  addition  to  the  measurements,  IDA  requires  initial  estimates  at 
both  the  object  and  PSF  as  well  as  support  regions  for  both  object  and  PSF  in 
order  to  begin  bhnd  deconvolution  processing.  As  alluded  to  in  Section  3.7.2,  three 
different  support  regions  are  selected  for  each  object  which  utihze  some  of  the  known 
attributes  of  the  true  impulse  response. 
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Assuming  circular  symmetry  in  the  PSF,  a  circle  is  chosen  as  the  shape  of  the 
support  region  for  the  PSF  estimate.  Knowing  the  objects  are  point  sources  which 
image  into  an  impulse  response  or  sum  of  impulse  responses,  support  regions  for 
the  objects  are  also  selected  to  be  circular  in  shape.  Therefore,  depending  upon  the 
number  of  stars  in  the  image,  the  support  region  for  the  object  estimate  consists  of  the 
total  area  under  separate  circular  regions  centered  over  each  “known”  star  location. 
Use  of  this  a  prion  information  tests  the  best  reconstruction  IDA  can  obtain.  For 
actual  imagery,  star  locations  may  be  revealed  from  other  more  complicated  image 
processing  techniques.  In  that  case,  support  regions  may  be  set  up  using  the  known 
data  and  subsequent  images  processed  via  this  bhnd  deconvolution  technique. 


The  size  of  the  circular  regions  must  be  determined.  An  appropriate  support 
size  might  also  be  related  to  the  impulse  response.  In  order  to  tie  the  support  size  to 
the  PSF,  setting  the  radius  of  the  support  circles  to  the  first,  second  and  third  zeros 
of  the  PSF  is  proposed.  Since  the  true  PSF  in  this  case  is  the  Airy  disk,  these  zeros 


correspond  to  the  radial  points  where  the  irradiance  for  an  Airy  function  is  zero.  Sir 
George  BiddeU  Airy  (1801-1892)  first  derived  an  equation  defining  the  irradiance  of 
a  point  source  imaged  through  a  circular  aperture  (10).  Airy’s  function  for  intensity. 


I{0),  appears  below 


m  =  7(0) 


2Ji{ka  sin  9) 
ka  sin  9 


where  9  represents  the  angular  distance  from  the  center,  k  is  the  wave  number,  and 
Ji(-)  is  a  Bessel  function  of  the  first  kind  of  order  one.  The  first  three  zeros  for  the 
Airy  function  occur  when  the  quantity  kasin9  is  3.83,  7.02  and  10.17,  respectively. 
Upon  integrating  the  Airy  function,  one  finds  that  84%  of  the  light  energy  is  present 
within  the  first  zero  corresponding  to  the  main  lobe  of  the  Airy  function.  91%  of 
the  fight  energy  remains  within  the  second  dark  ring  and  over  95%  lies  within  the 
third  zero  (10).  These  zeros  correspond  to  pixel  locations  for  the  true  PSF  used 
in  this  investigation.  Figure  4.1  shows  a  slice  through  the  center  of  the  true  PSF 
used  to  blur  the  objects  for  this  experiment.  Zeros  for  this  PSF  occur  at  11,  20  and 
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Point  Spread  Function  Slice 


Figure  4.1  A  slice  through  the  true  point  spread  function  reveals  zeros  at  11,  20 
and  29  pixels  from  the  center. 

29  pixels  from  the  center.  These  pixel  values  then  become  the  radii  of  the  support 
region  circles  used  to  create  support  arrays  for  the  object  and  PSF  estimates. 

In  this  first  experiment,  IDA  mimics  Lane’s  method  for  bhnd  deconvolution 
since  only  the  image  and  convolution  error  terms  contribute  to  the  combined  error 
function.  Additionally,  there  is  no  need  to  implement  Jefferies’  band-hmit  mask  as 
part  of  the  convolution  error  since  no  information  in  the  measured  image  spectrum 
exists  outside  of  the  cutoff  frequency  defined  by  the  OTF.  Thus,  applying  this  mask 
has  no  effect  on  the  noise-free  images  processed  in  Experiment  1.  A  total  of  twelve 
runs  of  the  deconvolution  algorithm  are  required  in  this  case  to  process  the  four 
images  for  each  of  three  different  support  regions. 

4-2.2  Simulation  Results.  Results  from  the  computer  simulation  runs  us¬ 
ing  unconstrained  minimization  for  noise-free  images  appear  outstanding.  Figure  4.2 
displays  a  gray-scaled  image  of  the  binary  star  with  the  smallest  separation — blurred 
Object  D,  now  called  Image  D.  After  bhnd  deconvolution  of  that  measurement  using 
the  Lane  method,  the  object  estimate  clearly  reveals  the  binary  star  as  shown  in  Fig- 
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Image  D 


Figure  4.2  Image  D,  a  noise-free  blurred  image. 


ure  4.3.  Three-dimensional  plots  of  the  object  and  PSF  estimates  resulting  from  the 
processed  images  appear  in  Appendix  C,  Figures  C.l  through  C.24  for  Experiment 
1. 

Section  3.6  defined  several  metrics  which  quantify  the  errors  in  the  object  and 
PSF  estimates  in  relation  to  the  actual  objects  and  PSF.  The  numerical  quantities 
appear  in  Table  4.1.  Each  of  the  four  images  were  processed  using  each  of  three 
separate  support  regions  defined  at  either  the  first,  second  or  third  zero  in  the  true 
PSF.  The  table  refers  to  the  different  support  regions  as  1,  2  or  3,  where  1,  2  and 
3  represent  a  circular  support  radius  equal  to  the  distance  to  the  first,  second  and 
third  zero  in  the  true  PSF,  respectively.  Also  hsted  are  the  number  of  iterations 
required  by  the  Lane  method  to  achieve  a  local  minimum.  The  total  errors  for  both 
the  object  and  the  PSF  (^;i)  represent  the  deviation  from  the  true  object  or  PSF 
within  the  passband  as  defined  by  Equations  3.6  and  3.7.  The  location  error  reveals 
the  ratio  of  the  error  in  the  point  source  location  to  the  blur  radius.  In  the  noise-free 


4-4 


Object  D  Estimate 


Figure  4.3  Object  estimate  of  Image  D  after  deconvolution  using  a  support  region 
extending  out  to  the  second  zero  of  the  true  PSF. 


case,  almost  all  point  sources  were  estimated  to  be  at  the  exact  same  pixel  location 
as  in  the  true  object,  thus  having  zero  location  error.  The  intensity  ratio  appears 
for  comparison  against  the  true  objects  B,  C  and  D — each  containing  a  binary  star 
with  an  intensity  ratio  of  1:0.7. 

For  further  analysis  of  the  support  constraint,  plots  are  shown  which  reveal  the 
spatial  frequency  content  of  the  true  objects  and  their  respective  estimates.  Figures 
4.4  -  4.7  present  the  radially  averaged,  normalized  object  spectra  estimated  using 
differing  support  regions.  In  each  plot,  the  horizontal  axis  displays  spatial  frequency 
where  “1”  represents  the  cutoff  frequency  of  the  optical  system.  The  vertical  axis  is 
normalized  for  comparison  purposes,  and  the  magnitudes  of  the  object  spectra  are 
radially  averaged  from  the  {u,v)  =  (0,0)  or  DC  value.  Graphed  are  the  following 
five  functions  within  each  figure. 

1.  First,  the  true  object  spectrum  appears.  This  spectrum  represents  the  single 
or  double  point  sources;  therefore,  it  has  infinite  frequency  content. 
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Table  4.1  Error  metric  summary  for  Experiment  1. 


OBJECT  /  Support 

Iter  “ 

0 

600' 

1:7/“= 

A/1 

294 

0.258 

0.279 

0.000 

n/a 

A  /  2 

564 

0.111 

0.060 

0.000 

n/a 

A  /  3 

873 

0.116 

0.045 

0.000 

n/a 

B/1 

3000 

1.155 

1.541 

0.545  /  0.545 

1:0.69 

B  /  2 

1184 

0.096 

0.046 

0.000  /  0.000 

1:0.70 

B  /  3 

1533 

0.146 

0.054 

0.000  /  0.000 

1:0.69 

C/1 

3000 

0.987 

1.867 

*  *  * 

C  /  2 

948 

0.097 

0.051 

0.000  /  0.000 

1:0.70 

C  /  3 

1057 

0.114 

0.041 

0.000  /  0.000 

1:0.71 

D/1 

1906 

0.661 

1.195 

0.455  /  0.530 

1:0.74 

D  /  2 

1140 

0.083 

0.051 

0.000  /  0.000 

1:0.72 

D  /  3 

1052 

0.102 

0.043 

0.000  /  0.000 

1:0.68 

“Iterations  required  by  IDA  for  reconstruction. 

^’Location  errors  are  listed  for  the  number  of  point  sources  present  in  the  object. 
“Intensity  ratio  is  only  applicable  for  binary  stars, 

“^For  this  estimate,  the  point  source  locations  and  relative  intensities  were  indiscernible. 


2.  The  solid  line  represents  the  true  object  spectrum  passed  through  a  lowpass 
filter  with  a  cutoff  frequency  equal  to  that  of  the  blurring  function.  This  soHd  hne 
reflects  the  theoretically  best  achievable  solution  possible,  given  that  all  information 
outside  of  the  band-hmit  was  discarded  by  the  imaging  process. 

3.  The  line  furthest  away  from  the  desired  spectrum,  shown  by  a  dotted  hne  in 
aU  four  graphs,  represents  the  estimates  created  using  the  smallest  support  region — 
defined  by  the  first  zero  in  the  PSF.  Inclusion  of  only  the  main  lobe  of  the  PSF 
in  the  support  region  does  not  provide  adequate  information  to  the  algorithm  to 
reconstruct  the  objects. 

4.  AU  four  estimates  generated  from  the  support  region  using  the  second  zero 
in  the  PSF  come  closest  to  the  low-pass  filtered  true  object  spectrum.  This  observed 
quahty  in  the  object  spectra  is  revealed  quantitatively  in  Table  4.1,  where  the  total 
object  error,  is  optimal  for  support  number  2  for  each  of  the  four  images. 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.4  Experiment  1;  Radially  averaged  object  spectra  for  Image  A  comparing 
the  filtered  true  object  spectrum  with  three  estimates  using  different 
support  regions. 

5.  Object  spectra  representing  the  estimates  generated  from  the  support  region 
using  the  third  zero  in  the  PSF  remain  close  to  the  estimates  processed  with  support 
region  2.  However,  the  reconstructions  shown  in  the  spatial  frequency  domain  reflect 
a  shghtly  diminished  tendency  toward  following  the  true  object  spectrum  within  the 
passband. 

4-2.3  Conclusions.  The  object  and  PSF  estimates  shown  in  Appendix  C, 
the  comparison  metrics  in  Table  4.1,  and  the  data  revealed  in  Figures  4.4  through  4.7 
indicate  that  Lane’s  unconstrained  minimization  approach  to  blind  deconvolution 
is  certainly  an  acceptable  technique  for  reconstructing  blurred  noise-free  images. 
Obviously,  an  accurate  estimate  of  the  support  for  both  the  object  an  PSF  aids  in 
this  reconstruction.  For  the  noise-free,  point  source  objects  apphed  in  Experiment 
1,  a  support  region  which  extends  out  to  the  second  zero  in  the  PSF  provides  the 
best  estimates.  This  “second  zero”  support  region  will  be  utiUzed  exclusively  for 
subsequent  experimentation.  The  high  frequency  content  in  the  spectra  and  “spikes” 
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Normalized  Ol^ed  Spectra  (radially  averaged) 


Figure  4.5  Experiment  1:  Radially  averaged  object  spectra  for  Image  B  comparing 
the  filtered  true  object  spectrum  with  three  estimates  using  different 
support  regions. 


Normatized  Object  Spectra  (radially  averaged) 


Figure  4.6  Experiment  1:  Radially  averaged  object  spectra  for  Image  C  comparing 
the  filtered  true  object  spectrum  with  three  estimates  using  different 
support  regions. 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.7  Experiment  1:  Radially  averaged  object  spectra  for  Image  D  comparing 
the  filtered  true  object  spectrum  with  three  estimates  using  different 
support  regions. 


visible  in  Figure  4.3  at  regular  intervals  as  well  as  the  other  object  estimates  deserve 
further  discussion. 

Both  the  high  spatial  frequencies  and  the  ring  of  “spikes”  probably  have  a 
great  deal  of  correlation,  especially  since  the  spikes  have  infinite  frequency  content. 
The  spikes  visible  in  all  of  the  object  estimates  shown  in  Appendix  C  are  formed  on 
the  very  edge  of  the  support  region,  thus  viewing  an  estimate  gives  one  a  feel  for 
the  size  of  the  support  region  applied.  The  boundary  of  the  support  region  presents 
a  large  discontinuity  in  the  image  plane.  This  discontinuity  may  be  the  single  most 
important  factor  in  understanding  the  existence  of  the  spikes.  Since  the  appearance 
of  the  spikes  occur  every  time,  one  might  reason  that  their  ehmination  by  filtering 
out  all  irradiance  present  on  the  edge  of  the  support  region — essentially,  erasing  the 
spikes — would  be  valid.  Another  viable  option  would  place  the  measured  image  in 
the  center  of  a  much  larger  array,  then  choose  a  support  region  large  enough  that 
the  discontinuity  occurs  outside  the  original  array  size.  This  avenue  is  alluded  to  in 
Jefferies’  statement  that  his  “support  constraints  were  very  loose”  (12). 
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Spatial  frequencies  that  are  estimated  above  the  cutoff  frequency  of  the  opti¬ 
cal  system  pose  more  questions.  The  selected  means  of  unconstrained  minimization 
imposes  no  strict  requirements  on  either  object  or  PSF  estimate.  Neither  the  im¬ 
age  error  nor  the  convolution  error  restrict  spatial  frequency  content  outside  the 
optical  band-limit;  thus,  Lane’s  method  will  always  produce  information  above  the 
band-hmit.  Jefferies  and  Christou  claim  that  super-resolution  is  possible  using  the 
unconstrained  minimization  technique.  Jefferies  cites  Lucy  [1992]  who  demonstrated 
that  if  the  convolution  components  are  of  finite  extent,  then  the  positivity  constraint 
permits  gains  in  resolution  over  the  diffraction  hmit  (12).  Goodman  [1968]  also  states 
that  retrieving  object  spectrum  values  outside  the  passband  is  theoretically  possi¬ 
ble  using  analytic  continuation  as  long  as  the  object  is  spatially  bounded  and  the 
spectral  information  within  the  passband  is  known  exactly  (7).  Others  refute  this 
super-resolution  claim,  stating  that  all  information  beyond  the  cutoff  frequency  of 
the  imaging  system  is  irrecoverable  (2).  Regardless  of  the  super-resolution  debate, 
high  frequency  content  present  in  the  object  spectrum  gives  rise  to  sharp  edges  within 
the  object  in  addition  to  the  spikes  present  around  the  support  edge.  Thus,  it  re¬ 
mains  to  be  proved  elsewhere  whether  the  discontinuity  at  the  edge  of  the  support 
region  creates  the  spikes  or  whether  the  high  frequency  content  allowed  to  exist  in 
the  object  spectrum  results  in  forming  the  spikes  around  the  support  region. 

4-3  Experiment  Two 

Determine  the  effect  of  photon  noise  on  the  object  estimate  and  compare  with 
results  previously  obtained  by  Lane.  Section  3.4  defines  the  objective  of  Experiment 
2  for  applying  the  unconstrained  minimization  technique  modified  by  Jefferies  to 
images  corrupted  by  photon  noise.  Jefferies’  modification  involves  the  addition  of 
a  band-hmit  filter  to  the  convolution  error  term  in  Equation  3.4  which  means  that 
the  convolution  of  the  object  and  PSF  estimate  do  not  include  the  high  frequency 
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noise  present  in  the  measured  image.  These  reconstructions  can  then  be  compared 
against  noisy  images  reconstructed  using  the  strict  Lane  method  without  the  mask. 

4.3.1  Simulation  Parameters.  The  only  variable  parameter  for  this  ex¬ 
periment  is  the  light  level  present  in  each  simulated  image.  Section  3.7.3  describes 
the  Poisson  model  used  to  generate  simulated  low  light  level  images.  Experiment  2 
utilizes  the  same  four  blurred  images  from  experiment  one  but  corrupts  each  with 
differing  amounts  of  photon  noise.  Photon  noise  is  quantified  in  terms  of  an  average 
number  of  photo  events  present  in  each  image.  For  this  experiment,  average  photo 
events  per  image,  K ,  of  100,000,  10,000,  and  5,000  are  used.  The  resulting  noisy 
images  are  shown  in  Appendix  B,  Figures  B.IO  through  B.21.  These  twelve  images 
become  the  measured  images  input  into  IDA  for  subsequent  processing. 

To  maintain  consistency  among  the  images,  a  single  support  region  correspond¬ 
ing  to  the  second  zero  in  the  true  PSF  is  utilized  throughout  this  experiment.  Ad¬ 
ditionally,  a  band-limit  equal  to  the  cutoff  frequency  in  the  OTF  used  to  blur  the 
objects  is  required  to  test  the  modification  to  the  Lane  method  recommended  by 
Jefferies  and  Christou.  Finally,  as  in  Experiment  1,  Experiment  2  only  applies  the 
image  and  convolution  error  to  the  combined  error  function  for  minimization. 

4.3.2  Simulation  Results.  Appendix  D  presents  the  actual  object  estimates 
obtained  in  Experiment  2  while  Table  4.2  summarizes  the  comparison  metrics  for 
the  object  and  PSF  estimates. 

The  first  two  rows  listed  for  each  separate  image  represent  a  comparison  be¬ 
tween  Lane’s  method  and  the  band-limit  modification.  Both  runs  used  the  image 
containing  10®  photo  events  as  the  measured  image.  The  values  for  the  object  error, 
and  the  PSF  error,  (h,  reveal  a  significant  decrease  in  the  error  value  when  the 
band-limit  mask  is  employed.  Further  evidence  can  be  observed  by  comparing  the 
resulting  object  estimates  visually  as  presented  in  Appendix  D,  Figures  D.l  through 
D.8.  Figures  D.l,  D.3,  D.5  and  D.7  show  the  object  estimates  using  the  strict  Lane 
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Table  4.2  Error  metric  summary  for  Experiment  2. 


OBJECT 

Egg 

mm 

6oc 

A 

■EKfl 

769 

0.437 

0.182 

n/a 

A 

875 

0.207 

0.203 

n/a 

A 

0.264 

0.288 

n/a 

A 

1347 

0.179 

0.152 

0.328 

n/a 

B 

0.091  /  0.091 

1:0.69 

B 

0.048 

0.000  /  0.091 

1:0.77 

B 

2743 

0.138 

0.077 

0.091  /  0.091 

1:0.76 

B 

1530 

0.292 

0.412 

0.288  /  0.257 

1:0.77 

C 

ww 

832 

0.316 

0.302 

0.257  1  *  *  * 

:(<  He  *  “ 

C 

1995 

0.080 

0.041 

0.091  /  0.128 

■Era 

C 

1109 

0.193 

0.155 

0.182  /  0.091 

C 

Imim 

1643 

0.313 

0.460 

0.364  /  0.203 

mi 

D 

mjjm 

1042 

0.258 

0.235 

^ 

D 

1270 

0.105 

0.066 

0.091  /  0.128 

lEia 

D 

1221 

0.292 

0.392 

0.182  /  0.203 

lira 

D 

3000 

0.595 

0.691 

0.386  /  0.273 

IB 

“Iterations  required  by  IDA  for  reconstruction. 

^Location  errors  are  listed  for  the  number  of  point  sources  present  in  the  object. 

“Intensity  ratio  is  only  applicable  for  binary  stars. 

“^The  first  line  for  each  image  [K  =  10®)  displays  the  metrics  for  estimates  resulting  from  the 
strict  Lane  method.  Subsequent  runs  utilize  the  band-limit  mask. 

“The  dimmer  point  source  location  was  indiscernible.  Thus,  no  intensity  ratio  was  calculated. 
■^See  footnote  e. 


method  on  the  photon  noise  corrupted  Images  A,  B,  C  and  D,  respectively.  Figures 
D.2,  D.4,  D.6  and  D.8  reveal  the  smooth  nature  of  the  object  estimates  utihzing  the 
band-hmit  mask.  The  Lane  method  proved  unsuccessful  at  reconstructing  images 
containing  fewer  than  10®  photo  events.  Therefore,  the  only  comparison  made  was 
for  images  where  .^  =  10®. 

Another  method  of  viewing  the  results  of  the  band-hmit  mask  comparison 
involves  plotting  the  radially  averaged  object  estimate  spectra.  Figures  4.8  through 
4.11  plot  these  spectra  along  with  the  filtered  and  unfiltered  true  object  spectra.  As 
previously  determined,  the  estimates  made  using  the  band-hmit  mask  more  closely 
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Normalized  Object  Spectra  (radially  averaged) 


Spatial  Frequency  (1  a  cutoff  freq) 

Figure  4.8  Experiment  2:  Radially  averaged  object  spectra  for  Image  A  [K  =  10®) 
comparing  the  use  of  the  band-hmit  mask  to  Lane’s  method.  The  filtered 
and  unfiltered  true  object  spectra  appear  as  reference. 

follow  the  true  object  spectrum  than  those  using  Lane’s  method.  Evidence  from  the 
spectra  as  well  as  the  object  estimate  lead  one  to  conclude  that  the  high  frequency 
noise  content  in  the  image  spectrum  hinders  the  reconstruction  process. 

Knowing  that  the  band-hmit  mask  on  the  convolution  error  aides  in  the  recon¬ 
struction  of  noisy  images,  results  from  images  with  differing  hght  levels  can  now  be 
presented.  Figures  D.IO  through  D.17  in  Appendix  D  present  the  object  estimates 
resulting  from  average  photo  event  levels  of  10,000  and  5,000.  As  previously  shown 
by  Lane  (13),  photon  level  significantly  affects  the  reconstruction  of  images  Hmited 
by  the  length  of  the  exposure  or  the  object’s  intensity.  Table  4.2  reveals  that  in  al¬ 
most  aU  cases,  decreasing  the  average  number  of  photo  events  increases  the  error  in 
both  the  object  and  PSF  estimates.  The  same  conclusion  is  reached  when  one  views 
the  radially  averaged  object  spectra  shown  in  Figures  4.12  through  4.15.  Except  for 
Image  A,  each  spectrum  steps  further  away  from  the  filtered  true  object  spectrum 
as  K  decreases.  As  the  number  of  photo  events  decreases,  the  noise  present  in  the 
measurement  increases  across  spectral  Hnes.  Therefore,  the  increased  noise  within 
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Nomnalized  Object  Spectra  (radially  averaged) 


Figure  4.9  Experiment  2:  Radially  averaged  object  spectra  for  Image  B  {K  =  10®) 
comparing  the  use  of  the  band-limit  mask  to  Lane’s  method.  The  filtered 
and  unfiltered  true  object  spectra  appear  as  reference. 


Normalized  Object  Spectra  (radially  averaged) 


Figure  4.10  Experiment  2;  Radially  averaged  object  spectra  for  Image  C  {K  =  10®) 
comparing  the  use  of  the  band-limit  mask  to  Lane’s  method.  The 
filtered  and  unfiltered  true  object  spectra  appear  as  reference. 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.11  Experiment  2:  Radially  averaged  object  spectra  for  Image  D  {K  =  10®) 
comparing  the  use  of  the  band-limit  mask  to  Lane’s  method.  The 
filtered  and  unfiltered  true  object  spectra  appear  as  reference. 


the  passband  manifests  itself  as  greater  error  within  the  object  and  PSF  estimates. 
For  images  with  around  1,000  photo  events,  the  resulting  object  and  PSF  estimates 
are  very  poor  and  thus,  not  presented  in  this  analysis. 


4.3.3  Conclusions.  Experiment  2  was  designed  to  test  the  validity  of  the 
band-limit  mask  as  applied  to  photon  limited  images  as  well  as  determine  the  ability 
of  IDA  to  reconstruct  images  with  low  light  levels.  As  expected,  decreasing  the 
light  level  or  number  of  photo  events  present  in  an  image  significantly  degrades  the 
reconstruction.  The  random  arrival  of  photo  events  within  an  image  effectively  act 
as  noise.  Therefore,  a  limit  to  the  reconstruction  is  reached  when  the  photon  noise 
sufficiently  conceals  the  signal  present  in  the  image. 

The  results  definitively  prove  the  positive  utility  of  the  band-limit  mask  on 
the  convolution  error.  This  mask  allows  Jefferies’  algorithm  to  reconstruct  images 
with  far  less  light  present  in  the  measurement  than  with  Lane’s  method.  The  mask 
effectively  filters  out  the  high  frequency  content  of  the  noise,  and  subsequently  allows 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.12  Experiment  2:  Radially  averaged  object  spectra  for  Image  A  applying 
the  band-Hmit  mask  to  images  with  differing  amounts  of  photon  noise, 
K  =  10®,  10“*  and  5, 000.  The  filtered  and  unfiltered  true  object  spectra 
appear  as  reference. 


Normalized  Object  Spectra  (radially  averaged) 


Figure  4.13  Experiment  2:  Radially  averaged  object  spectra  for  Image  B  applying 
the  band-hmit  mask  to  images  with  differing  amounts  of  photon  noise, 
K  =  10®,  10^  and  5, 000.  The  filtered  and  unfiltered  true  object  spectra 
appear  as  reference. 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.14  Experiment  2:  Radially  averaged  object  spectra  for  Image  C  applying 
the  band-limit  mask  to  images  with  differing  amounts  of  photon  noise, 
K  =  10®,  10"^  and  5, 000.  The  filtered  and  unfiltered  true  object  spectra 
appear  as  reference. 


Normalized  Object  Spectra  (radially  averaged) 


Figure  4.15  Experiment  2:  Radially  averaged  object  spectra  for  Image  D  applying 
the  band-limit  mask  to  images  with  differing  amounts  of  photon  noise, 
=  10®,  10“*  and  5, 000.  The  filtered  and  unfiltered  true  object  spectra 
appear  as  reference. 
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reconstruction  of  images  with  a  much  lower  signal-to-noise  ratio.  As  revealed  in  Lane 
[1992],  unconstrained  minimization  cannot  effectively  reconstruct  images  when  the 
number  of  photo  events  in  the  brightest  pixel  is  less  than  10“*  (13).  With  the  band- 
Hmit  mask,  however,  excellent  results  can  be  obtained  when  the  brightest  pixel  has 
less  than  10^  photo  events  present  as  expressed  in  Figure  D.14  which  shows  the 
reconstruction  of  Figure  B.17  containing  fewer  than  70  photo  events  in  its  brightest 
pixel. 


4-4  Experiment  Three 

Determine  the  utility  of  the  band-limit  error  term  toward  improving  the  recon¬ 
struction  process.  With  the  success  of  the  band-limit  mask  in  the  convolution  error 
term,  this  final  experiment  studies  the  band-hmit  error  term,  shown  in  Equation 
3.5,  which  requires  a  knowledge  of  the  same  information  necessary  in  Experiment  2. 
If  the  cutoff  frequency  of  the  optical  system  is  known  for  a  given  image,  then  the 
band-limit  error  term  might  prove  useful  in  the  reconstruction  process  just  as  the 
band-hmit  mask  greatly  improved  the  object  estimates  in  the  previous  experiment. 
A  look  at  the  comparison  metrics  and  the  number  of  iterations  required  in  this  ex¬ 
periment  compared  with  similar  quantities  available  from  previous  experiments  for 
the  same  measurement  should  provide  an  understanding  of  the  effect  the  band-hmit 
error  has  on  the  bhnd  deconvolution  process. 

4-4-i  Simulation  Parameters.  Experiment  3  differs  from  previous  runs 
using  this  unconstrained  minimization  technique  in  that  a  third  error  term  is  added 
to  the  objective  function  to  be  minimized.  The  only  input  this  third  term  requires  is 
the  band-hmit  of  the  imaging  system,  which  is  also  required  for  use  of  the  band-hmit 
mask  in  the  convolution  error  term.  The  convolution  error  term  retains  the  band- 
hmit  mask  and  thus  remains  unchanged  from  Experiment  2.  The  image  domain 
error  term  stiU  minimizes  any  negative  values  in  either  the  object  or  PSF  estimates. 
As  for  support,  the  optimal  support  size  determined  in  Experiment  1  wiU  be  input 
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once  again  in  this  study.  A  priori^  one  must  know  the  linear  extent  of  the  second 
zero  in  the  impulse  response  in  the  image  domain  and  know  the  band-Hmit  of  the 
impulse  response  in  the  frequency  domain. 

Selected  images  used  in  this  experiment  include  the  four  infinite  signal  blurred 
images  applied  in  Experiment  1  and  their  simulated,  light  limited  images  used  in 
Experiment  2  with  average  photo  events  of  10®  and  10^.  This  experiment  wiU  input 
twelve  different  simulated  images  into  the  iterative  deconvolution  algorithm  to  test 
the  utility  of  the  band-Hmit  error  term.  Jefferies  and  Christou  state  that  the  use 
of  the  band-limit  term  must  “be  relaxed  before  convergence”  (12).  Therefore,  one 
can  expect  that  the  algorithm  wiU  run  for  a  number  of  iterations  with  the  band- 
Hmit  error  term,  then  converge  after  a  series  of  more  iterations  using  only  the  image 
domain  and  convolution  error  terms. 

Some  preparatory  analysis  provides  useful  data  regarding  the  use  of  the  band- 
Hmit  term.  The  purpose  for  this  term  is  to  reduce  high  spatial  frequency  content  in 
the  PSF  estimate.  Figures  4.16  and  4.17  reveal  the  consequence  of  allowing  the  band- 
Hmit  error  to  be  used  until  convergence  for  noise-free  Image  B.  The  PSF  estimate  is 
very  smooth  and  has  the  appearance  of  an  Airy  function,  yet  its  error,  ^/i,  of  0.423  is 
much  larger  than  the  error  value  obtained  in  Experiment  1  for  the  same  measurement 
and  support — see  Image  B,  support  2  in  Table  4.1.  On  the  other  hand,  the  object 
estimate  consists  of  two  rings  centered  on  the  actual  point  source  locations.  The 
object  error,  is  0.376  which  is  almost  four  times  greater  than  the  error  obtained 
in  Experiment  1  (see  Table  4.1).  As  Figure  4.17  shows,  only  the  rings  contain  the 
intensity  and  the  points  inside  the  ring  have  no  intensity.  Obviously,  allowing  the 
band-Hmit  term  to  remain  in  the  error  function  until  convergence  increases  the  object 
and  PSF  errors.  The  method  used  for  this  experiment  allows  the  band-Hmit  term  to 
be  used  for  the  first  fifty  iterations,  then  the  current  estimates  at  iteration  fifty  are 
input  as  the  initial  guesses  for  the  algorithm  using  only  the  image  domain  and  the 
convolution  error  terms. 
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PSF  estimate 


Figure  4.16  Experiment  3:  PSF  estimate  of  Image  B  after  2693  iterations  using  the 
band-limit  error  term  to  convergence.  The  band-limit  term  effectively 
constrains  high  spatial  frequencies  within  the  PSF  estimates. 


Object  estimate 


Figure  4.17  Experiment  3;  Object  estimate  of  Image  B  after  2693  iterations  using 
the  band-limit  error  term  to  convergence.  Note  the  circular  rings  which 
surround  the  true  point  source  locations,  though  these  locations  possess 
no  intensity. 
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4.4- ^  Simulation  Results.  Strict  application  of  the  band-limit  error  gives 
inconsistent  results.  Comparing  these  results  with  those  estabhshed  in  previous 
experiments  reveal  a  wide  range  of  differences.  Figures  4.18  through  4.21  show  the 
results  of  processing  the  K  —  10®  images  with  the  band-limit  error  next  to  the  results 
for  the  same  simulated  image  in  Experiment  2.  Again,  the  filtered  and  true  object 
spectra  appear  as  references.  From  the  graph,  Image  D  (Figure  4.21)  shows  httle 
deviation  in  the  two  reconstructions  and  would  lead  one  to  beheve  that  the  band- 
hmit  error  has  minimal  effect  on  the  algorithm.  However,  Image  A  (Figure  4.18) 
reveals  an  improvement  with  the  band-limit  error  whereas  Images  B  and  C  (Figures 
4.19  and  4.20)  show  that  appUcation  of  the  error  has  a  detrimental  effect  on  the 
object  estimates. 

Table  4.3  presents  the  comparison  metrics  for  each  object  and  PSF  estimate. 
Metrics  for  the  infinite  signal  images  can  be  compared  with  data  from  Experiment  1 
using  the  same  support  region.  The  use  of  the  band-limit  error  term  improves  Image 
A,  but  requires  almost  four  times  as  many  iterations.  Image  D  has  very  similar 
results  in  both  cases,  but  metrics  for  Images  B  and  C  were  much  lower  without  the 
band-Hmit  error  term  in  Experiment  1  as  shown  in  Table  4.1. 

Similar  results  were  obtained  applying  images  corrupted  by  photon  noise  to 
this  study.  When  compared  against  similar  data  in  Table  4.2  for  the  simulated 
photon  Hmited  measurements  from  Experiment  2,  the  band-limit  error  significantly 
improves  Image  A  in  terms  of  the  metrics,  and  ^h.  However,  this  error  shows 
little  deviation  from  previous  results  for  Images  C  and  D  and  proves  detrimental  for 
the  reconstruction  of  Image  B.  Appendix  E,  Figures  E.l  through  E.12,  displays  the 
resulting  object  estimates  for  all  of  the  data  presented  in  Table  4.3. 

4.4- 3  Conclusions.  The  inconsistent  results  make  drawing  conclusions  from 
Experiment  3  very  difficult.  From  the  limited  number  of  images  processed,  the  best 
results  in  using  the  band-limit  error  occur  when  the  point  sources  are  grouped  closely 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.18  Experiment  3:  Radially  averaged  object  spectra  for  Image  A  comparing 
the  object  estimate  reconstructed  with  the  aide  of  the  band-limit  error 
to  the  estimate  found  in  experiment  2  without  the  error. 


Normalized  Object  Spectra  (radially  averaged) 


Figure  4.19  Experiment  3:  Radially  averaged  object  spectra  for  Image  B  comparing 
the  object  estimate  reconstructed  with  the  aide  of  the  band-limit  error 
to  the  estimate  found  in  experiment  2  without  the  error. 
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Normalized  Object  Spectra  (radially  averaged) 


Figure  4.20  Experiment  3:  Radially  averaged  object  spectra  for  Image  C  comparing 
the  object  estimate  reconstructed  with  the  aide  of  the  band-hmit  error 
to  the  estimate  found  in  experiment  2  without  the  error. 


Normalized  Object  Spectra  (radially  averaged) 


Figure  4.21  Experiment  3:  Radially  averaged  object  spectra  for  Image  D  comparing 
the  object  estimate  reconstructed  with  the  aide  of  the  band-hmit  error 
to  the  estimate  found  in  experiment  2  without  the  error. 
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Table  4.3  Error  metric  summary  for  Experiment  3. 


OBJECT 

K 

Iter  “ 

0 

mm 

Cloc  ^ 

1:7/“ 

A 

(X 

likn^ 

n/a 

A 

n/a 

A 

2739 

0.071 

n/a 

B 

CXD 

2693 

IBB 

0.635 

J)C  *  *  j  *  * 

*  *  *  ® 

B 

10® 

1625 

0.678 

0.328  /  0.328 

mm 

B 

10^* 

2169 

2.262 

0.257  /  0.257 

ISli 

C 

oo 

2113 

0.508 

He  >|c  itc  j  *  *  * 

*  *  s(t 

C 

10® 

2212 

0.091 

0.046 

m 

C 

10'* 

1569 

0.211 

0.185 

D 

oo 

0.083 

0.047 

0.091  /  0.128 

D 

10® 

IBSI 

0.063 

0.091  /  0.128 

D 

10'* 

0.415 

0.182  /  0.128 

“Iterations  include  fifty  with  the  band-limit  error  then  all  subsequent  iterations  without  the 
band-limit  error. 

^’Location  errors  are  listed  for  the  number  of  point  sources  present  in  the  object. 

“^Intensity  ratio  is  only  applicable  for  binary  stars. 

“^The  first  line  for  each  image  {K  =  oo)  displays  the  metrics  for  estimates  resulting  from  a 
noise-free,  infinite  signal  image. 

®The  dimmer  point  source  location  was  indiscernible.  No  intensity  ratio  was  calculated. 

■^See  note  e. 
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around  the  center  of  the  image.  For  Image  B  which  contains  two  point  sources  spaced 
far  apart,  unsatisfactory  results  were  obtained.  As  the  point  sources  approach  one 
another,  the  results  approach  those  of  Experiment  2.  These  trends  seem  to  indicate 
that  the  band-limit  error  works  better  when  applied  to  images  containing  information 
tightly  grouped  around  the  central  point  in  the  image  plane. 

In  terms  of  the  unconstrained  minimization  method,  the  addition  of  the  band- 
bmit  error  alters  the  objective  function  and  thus  changes  the  direction  of  the  search 
for  a  local  minimum.  In  some  cases,  this  change  in  search  direction  can  be  a  help 
and  in  others  it  might  be  a  hindrance.  Also  a  factor  may  be  the  number  of  iterations 
one  should  optimally  take  in  this  new  search  direction.  More  analysis  is  required 
to  determine  whether  the  number  of  iterations  performed  using  the  band-limit  error 
might  have  a  correlation  with  improved  results.  Perhaps  improvement  might  be  made 
for  different  images  using  a  different  number  of  initial  iterations  with  the  band-limit 
error. 

^.5  Summary 

This  investigation  sought  to  test  whether  unconstrained  minimization  might 
be  an  effective  tool  in  solving  the  blind  deconvolution  problem  for  compensated 
imagery.  Therefore,  experiments  were  run  on  simulated  blurred  images  with  and 
without  photon  noise  to  determine  the  abiUty  of  the  Lane  method  as  modified  by 
Jefferies  and  Christou  to  reconstruct  imagery.  Experimental  objectives  included  the 
following:  (1)  a  test  of  the  support  criterion,  (2)  analysis  of  the  band-limit  mask  in 
the  presence  of  photon  noise,  and  (3)  a  study  of  the  band-limit  error  term  used  by 
Jefferies  and  Christou.  Results  from  these  experiments  provided  answers  to  some 
questions,  but  generated  several  more.  For  instance,  the  high  spatial  frequency 
content  produced  in  the  noise-free  case  remains  in  the  object  and  PSF  estimates 
throughout  all  of  the  experiments.  Though  the  high  frequency  information  yields 
weU-defined  point  sources  apart  from  the  image  plane,  spikes  also  appear  on  the 
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edge  of  the  support  region.  The  vaHdity  of  the  high  frequency  information  also 
deserves  further  attention.  Analysis  may  prove  that  this  unconstrained  minimization 
technique  performs  some  analytic  continuation  as  described  by  Goodman  [1968]  as 
necessary  for  super-resolution  (7). 

Along  with  the  many  questions  evolving  from  this  investigation,  some  positive 
results  arose  also.  Resulting  data  proves  that  the  unconstrained  minimization  tech¬ 
nique  can  accurately  recover  object  information  from  a  measured  image  with  very 
Httle  information  regarding  the  impulse  response.  The  accuracy  of  the  object  and 
PSF  estimates  have  a  high  dependence  upon  the  choice  of  the  support  region.  Excel¬ 
lent  results  were  obtained  using  a  support  region  extending  out  to  the  second  zero  in 
the  true  point  spread  function.  Additionally,  the  band-hmit  mask  contained  in  the 
convolution  error  term  allows  the  algorithm  to  handle  more  noise  in  a  measurement 
than  the  strict  Lane  method.  This  mask  ensures  that  the  object  and  PSF  estimates 
convolve  to  form  the  portion  of  the  measurement  contained  within  the  passband  of 
the  optics.  Thus,  a  priori  information  regarding  the  cutoff  frequency  is  essential  in 
reconstructing  low  hght  level  images.  Lastly,  the  band-hmit  error  term  acts  as  an 
additional  factor  which  may  improve  the  quality  of  the  reconstructed  image.  Its  use 
win  prevent  the  trivial  solution  from  occurring,  but  using  the  band-hmit  term  to 
alter  the  search  pattern  of  the  minimization  routine  may  also  produce  some  positive 
effects  on  the  eventual  “best”  object  estimate. 
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V.  Conclusion  and  Recommendations 

5. 1  Introduction 

Imagery  obtained  through  the  adaptive  optics  telescope  at  AMOS  require  addi¬ 
tional  processing  to  remove  the  blurring  still  present  in  the  images.  Little  information 
is  known  about  the  blurring  function  which  produces  the  images  output  from  the 
adaptive  optics  telescope;  thus  a  means  of  processing  the  compensated  imagery  is 
necessary  which  does  not  require  knowledge  of  the  PSF.  The  term  given  to  this  type 
of  problem  is  blind  deconvolution. 

Over  the  last  few  years,  several  techniques  have  been  suggested  which  attempt 
to  iteratively  solve  the  bHnd  deconvolution  problem  with  subsequent  estimates  at  an 
object  and  a  PSF  under  certain  known  constraints.  Ayers  and  Dainty  developed  the 
first  working  algorithm  which  successively  constrained  object  and  PSF  estimates  in 
both  the  image  domain  and  spatial  frequency  domain.  The  Ayers-Dainty  technique 
also  had  some  problems  in  that  no  avenue  existed  within  the  algorithm  for  determin¬ 
ing  when  the  best  estimate  of  the  object  and  PSF  had  been  reached.  Holmes  and 
Lane  reached  separate  solutions  for  determining  convergence  criteria  for  the  bhnd 
deconvolution  problem.  Holmes  formulated  the  problem  in  terms  of  a  maximum- 
hkehhood  approach  and  Lane  developed  an  unconstrained  minimization  method. 
This  investigation  focused  on  analyzing  the  unconstrained  minimization  technique 
to  solve  the  problem  of  deconvolving  a  simulated  image  with  little  or  no  information 
about  the  impulse  response. 

Methods  were  explained  and  experiments  were  developed  around  analyzing  the 
Lane  method  to  solving  blind  deconvolution  as  modified  by  Jefferies  and  Christou 
[1993].  Simulated  imagery  was  used  to  test  the  response  of  the  iterative  deconvolu¬ 
tion  algorithm  to  noise-free  images  and  images  corrupted  exclusively  by  photon  noise. 
The  objectives  tested  involved  observing  the  results  of  differing  support  sizes,  study¬ 
ing  the  effects  of  a  band-limit  mask  to  reconstruct  noisy  images,  and  determining 
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the  utility  of  a  new  error  term  designed  to  reduce  high  spatial  frequency  content  in 
the  optical  transfer  function  estimates.  Results  produced  by  the  experiments  proved 
extremely  promising.  With  information  regarding  only  the  band-limit  of  the  optics, 
this  algorithm  effectively  deconvolved  noise-free  data  as  well  as  simulated  images 
modeled  to  contain  differing  amounts  of  photon  noise.  Several  important  conclu¬ 
sions  were  also  drawn  from  the  experimentation  that  deserve  further  enumeration. 

5.2  Conclusions 

In  addition  to  accurately  recovering  object  information  from  measurements 
with  an  unknown  impulse  response,  results  from  the  investigation  of  the  uncon¬ 
strained  minimization  technique  revealed  several  significant  conclusions. 

1.  The  accuracy  of  the  object  and  PSF  estimates  are  highly  dependent  upon 
the  choice  of  the  support  region.  Excellent  results  were  obtained  using  a  support 
region  extending  out  to  the  second  zero  in  the  PSF,  where  the  area  within  the  second 
zero  retains  91  percent  of  the  volume  of  the  actual  PSF. 

2.  The  band-limit  mask  contained  in  the  convolution  error  term  allows  the 
algorithm  to  handle  significantly  more  noise  in  a  measurement  than  the  strict  Lane 
method.  This  mask  ensures  that  the  object  and  PSF  estimates  convolve  to  form  the 
portion  of  the  measurement  contained  within  the  passband  of  the  optics.  Thus,  a 
priori  information  regarding  the  cutoff  frequency  is  essential  in  reconstructing  low 
light-level  images. 

3.  The  band-limit  error  term  acts  as  an  additional  factor  which  may  improve 
the  quality  of  reconstructed  images.  Its  use  wiU  always  prevent  the  trivial  solution 
from  occurring,  but  using  the  band-hmit  term  to  alter  the  search  pattern  of  the 
minimization  routine  may  also  produce  some  positive  effects  in  the  object  estimate. 
Since  the  band-limit  error  term  can  have  a  negative  impact  on  the  reconstruction, 
special  care  must  be  taken  with  regard  to  using  this  term. 
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4.  High  spatial  frequency  information,  outside  of  the  cutoff  frequency,  appears 
in  the  object  and  PSF  estimates.  Due  to  the  nature  of  unconstrained  minimization, 
nothing  prevents  the  inclusion  of  this  information.  The  information  manifests  itself 
in  the  object  estimate  as  well-defined  edges,  meaning  that  no  gradual  roU-off  occurs 
in  the  transition  from  high  intensity  regions  to  low  intensity  regions  within  the 
image  plane.  This  quality  improves  the  resolution  of  the  object  estimate  beyond 
the  diffraction  limit.  However,  the  high  frequencies  also  contribute  to  spikes  which 
form  at  regular  intervals  around  the  edge  of  the  support  region.  The  presence  of  this 
information  above  the  cutoff  frequency  may  lead  some  to  question  the  validity  of  the 
results. 

5.3  Recommendations  for  Further  Research 

Obviously  much  more  research  is  required  before  the  unconstrained  minimiza¬ 
tion  approach  to  blind  deconvolution  could  be  used  operationally  at  AMOS.  This 
investigation,  however,  seeks  to  present  a  basic  framework  to  be  built  upon  in  subse¬ 
quent  research  efforts.  A  single  blind  deconvolution  technique  is  analyzed  for  objects 
containing  single  and  double  point  sources.  The  following  items  are  recommenda¬ 
tions  for  future  research  aimed  at  producing  an  operational  algorithm  to  improve 
the  quality  of  compensated  imagery  similar  to  that  produced  at  AMOS. 

1.  Try  weighting  the  different  error  terms  in  the  objective  function.  Evidence 
from  Experiment  3  indicates  that  the  band-limit  error  term  can  help  or  hinder  the 
reconstruction  process.  A  variable  weight  on  the  band-limit  term  as  well  as  the 
image  domain  and  convolution  error  terms  may  reveal  an  improved  convergence  or 
possibly  enhanced  object  and  PSF  estimates. 

2.  Determine  the  steps  necessary  to  reconstruct  extended  objects  through  the 
unconstrained  minimization  technique.  Continuous  extended  objects,  like  satellites, 
pose  entirely  new  problems  when  imaging  through  the  atmosphere.  Although  out¬ 
standing  results  were  obtained  on  point  source  objects,  a  blurred  extended  object 
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may  prove  very  difficult  to  reconstruct  using  this  method.  The  process  of  testing 
the  iterative  deconvolution  algorithm  should  start  where  this  investigation  leaves 
off — with  multiple  point  source  objects.  Adding  more  point  sources  should  provide 
an  understanding  of  how  the  algorithm  performs  for  larger,  more  complex  images. 
Modifications  to  the  unconstrained  minimization  technique  are  most  hkely  required 
to  enable  the  recovery  of  detailed  information  in  extended  objects  through  bhnd 
deconvolution. 

3.  To  reconstruct  compensated  images  of  sateUites,  compare  unconstrained 
minimization  to  maximum-hkehhood  estimation  as  a  different  means  to  accomphsh 
the  goal  of  bhndly  deconvolving  extended  object  data.  A  more  thorough  analy¬ 
sis  of  the  work  by  Holmes  (11)  and  Schulz  (21)  may  lead  to  a  more  robust  blind 
deconvolution  algorithm  capable  of  handhng  the  complexity  of  extended  objects. 

4.  Analyze  the  source  of  the  high  spatial  frequency  information  found  in  the 
object  and  PSF  estimates  produced  by  this  algorithm.  Further  research  is  required 
on  Jefferies  and  Christou’s  claim  that  their  iterative  deconvolution  algorithm  can 
super-resolve  blurred  images  (12).  Such  research  would  certainly  help  validate  this 
method  and  give  credence  to  the  results  of  the  algorithm  which  contain  spatial 
frequency  information  beyond  the  band-limit  of  the  imaging  system. 


5-4 


Appendix  A.  Iterative  Deconvolution  Algorithm  Tutorial 

A .  1  Introduction 

This  investigation  focuses  on  the  iterative  deconvolution  algorithm  (IDA)  de¬ 
veloped  by  Stuart  Jefferies  and  Julian  Christou  who  graciously  provided  the  actual 
code  used  in  their  paper  (12).  The  IDA  code  (annotated  simply  IDA  for  this  ap¬ 
pendix)  is  an  extensive  image  processing  package.  It  provides  numerous  options  for 
performing  bUnd  deconvolution  and  “not-so-bhnd”  deconvolution.  The  program  can 
take  simply  a  measurement  and  initial  guesses  at  both  an  object  and  PSF  and  it¬ 
eratively  alter  the  object  and  PSF  estimates  toward  a  best  estimate  of  each  using 
the  Lane  method  (13)  of  unconstrained  minimization.  If  more  information  is  avail¬ 
able,  hke  the  band-hmit  of  the  optics  or  Fourier  modulus  data,  then  IDA  uses  this 
information  in  an  unconstrained  minimization  approach  which  is  a  modification  to 
the  Lane  method.  Additionally,  IDA  handles  multiple  images  of  the  same  object 
with  ease.  The  program  runs  on  a  SUN  workstation.  The  following  sections  discuss 
how  IDA  implements  the  algorithm  derived  by  Lane  and  modified  by  Jefferies  and 
Christou,  the  input  required  to  run  the  program  and  the  proper  format  for  the  data, 
what  options  are  available  in  processing  images,  and  the  modifications  which  must 
be  made  each  time  IDA  is  run  on  a  different  image. 

A. 2  How  IDA  works 

As  revealed  in  Section  3.5,  IDA  solves  the  blind  deconvolution  problem  through 
the  minimization  of  an  error  function.  Minimization  is  achieved  through  a  conjugate 
gradient  minimization  routine  which  finds  a  local  minimum  to  a  multi-dimensional 
objective  function  through  a  steepest-descent  search  method.  Such  a  method  requires 
only  a  subroutine  to  calculate  the  value  of  a  function  at  a  particular  location  and 
another  subroutine  to  calculate  the  derivative  at  the  location.  IDA  provides  these 
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subroutines  to  the  conjugate  gradient  minimization  function  using  the  error  terms 
and  their  derivatives  presented  in  Jefferies  and  Christou  [1993]. 

IDA  minimizes  the  combined  error  function  shown  in  Equation  2.16.  Each 
pixel  location  in  both  the  object  and  PSF  estimate  represent  the  different  variables 
from  which  the  error  function  must  be  minimized.  The  support  region  input  is 
utilized  by  the  program  to  reduce  the  number  of  variables.  IDA  only  includes  those 
pixel  locations  for  both  the  object  and  PSF  which  He  within  the  support  region.  AU 
other  pixels  are  set  to  zero  outside  the  support.  Using  the  arrays  provided  for  initial 
estimates,  IDA  picks  off  the  value  at  each  pixel  location  within  the  support  to  create 
the  series  of  values  for  a  first  guess  at  the  object  and  PSF  estimate. 

The  conjugate  gradient  minimization  routine  alters  each  pixel  value  in  the 
series  in  a  direction  opposite  to  the  gradient.  Thus,  each  successive  iteration  reduces 
the  combined  error  function  value.  The  iterations  cease  when  one  of  three  things 
occur:  (1)  a  local  minimum  is  reached  such  that  any  variation  in  pixel  values  causes 
an  increase  in  the  objective  function,  (2)  the  reduction  in  the  error  function  value 
for  each  iteration  decreases  below  a  pre-set  tolerance,  or  (3)  a  maximum  number  of 
iterations  is  reached. 

A. 3  Required  Input 

Running  IDA  requires  several  input  data  files.  Modifications  made  to  Jefferies 
and  Christou’s  code  for  this  investigation  require  that  the  data  files  be  formatted  in 
either  binary  or  ASCII.  Below,  each  of  the  inputs  are  defined  and  the  appropriate 
formats  are  annotated. 

A. 3.1  Convolution  Image.  IDA  first  requires  a  filename  containing  the 
data  representing  the  convolution  image,  also  known  as  the  measurement.  A  binary 
format  is  required  for  this  data. 
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A. 3. 2  Object  Estimate.  The  second  input  necessary  to  reconstruct  an 
object  from  the  input  measurement  is  the  filename  containing  the  initial  estimate 
at  the  object  in  binary  format.  At  regular  intervals  in  the  iteration  process,  the 
current  object  estimate  will  be  written  to  this  file.  A  binary  representation  of  the 
final  estimate  appears  in  this  file  alone  after  IDA  completes;  thus,  the  initial  estimate 
contained  in  this  file  is  destroyed. 

A. 3. 3  PSF  Estimate.  Similar  to  the  object  estimate,  IDA  requires  the  name 
of  the  binary  file  containing  the  first  guess  at  a  PSF.  This  estimate  is  also  overwritten 
by  successive  estimates.  Upon  completion  of  the  program,  this  file  contains  the  final 
estimate  at  the  PSF. 

A. 3. 4  Object  and  PSF  Support  Arrays.  A  single  file  containing  two  support 
arrays  must  be  provided  to  IDA.  The  file  contains  the  object  support  array  followed 
immediately  by  the  PSF  support  array  in  ASCII  format.  The  two  arrays  must  be 
the  same  size  as  the  object,  PSF  and  convolution  image.  The  value  in  each  pixel 
location  contains  either  a  “0”  or  “1”,  where  one  represents  a  pixel  within  the  support 
region  and  zero  represents  a  pixel  outside  the  support  region. 

A. 3. 5  Output  Log  File.  IDA  requires  a  unique  filename  to  write  the  output 
information.  A  name  of  an  existing  file  wiU  be  rejected.  This  file  wiU  contain  a  sum¬ 
mary  of  the  parameters  entered  to  run  the  program  and  a  hsting  of  the  quantitative 
values  of  the  different  error  terms  at  each  iteration. 

A. 3. 6  Maximum  Iterations.  The  final  input  required  for  IDA  to  reconstruct 
an  image  is  the  maximum  number  of  iterations  allowed  for  the  run.  Choosing  too 
few  iterations  is  not  a  problem  since  the  program  has  an  option  which  allows  the 
user  to  “restart”  from  the  point  at  which  the  maximum  number  of  iterations  was 
reached  by  inputting  the  estimates  obtained  as  the  initial  estimates  and  typing  “1” 
for  the  restart  parameter. 
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A. 4  Available  Options 

After  the  information  above  has  been  entered,  a  series  of  additional  questions 
regarding  which  error  terms  IDA  will  use  to  perform  unconstrained  minimization 
must  be  answered.  The  first  question  asks  for  either  0,  1,  or  2  for  the  telescope 
aperture  parameter.  The  second  queries  whether  a  restart  is  required,  and  the  third 
asks  for  a  Fourier  spectrum  parameter  of  0,  1,  or  2.  If  a  non-circular  convolution 
mask  must  be  input,  then  a  fourth  question  prompts  the  user  for  that  informa¬ 
tion.  The  next  few  sections  outline  how  to  implement  several  different  variations  on 
the  iterative  deconvolution  algorithm  by  answering  the  four  questions  listed  below 
differently. 

Ql.  Input  telescope  aperture  parameter  (2,1,0): 

Q2.  Input  1  for  a  restart  (else  0): 

Q3.  Input  Fourier  spectrum  parameter  (2,1,0): 

Q4.  Input  1  for  non-circular  convolution  mask  (else  0): 

A. 4-1  Basic  Lane  Method.  The  Lane  method  is  the  simplest  form  of  the 
unconstrained  minimization  technique  which  IDA  can  perform.  The  Lane  method, 
as  discussed  in  Section  2.5.3,  minimizes  an  error  function  containing  only  an  image 
domain  error  term  and  a  convolution  error  term.  Since  no  a  priori  information  is 
required,  a  “0”  is  entered  for  all  four  questions.  However,  if  a  restart  is  desired,  then 
“1”  would  be  the  appropriate  entry  for  the  second  question,  Q2. 

A. 4 -2  Band-limit  Mask.  As  shown  in  experiment  two,  the  addition  of  a 
band-hmiting  mask  on  the  convolution  error  term  significantly  increases  the  abihty 
of  IDA  to  reconstruct  images  in  the  presence  of  noise.  Utilization  of  the  band- 
limit  error  occurs  whenever  the  telescope  aperture  parameter,  set  in  Ql,  is  non-zero. 
Entering  either  “1”  or  “2”  for  Ql  will  require  an  additional  data  file.  After  all  the 
questions  are  asked,  IDA  will  prompt  the  user  for  the  file  containing  cutoff  frequency 
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information.  If  multiple  measurements  images  of  the  same  object  are  input,  then 
IDA  expects  to  find  the  cutoff  frequency  in  pixels  from  the  center  for  each  unknown 
PSF.  The  program  then  creates  a  binary  mask  array  for  each  PSF  with  a  circular 
region  of  “I’s”  with  a  radius  equal  to  that  input  from  the  cutoff  frequency  file.  If  non¬ 
circular  masks  are  desired,  then  answering  “1”  to  Q4  gives  the  user  the  opportunity 
to  enter  a  filename  containing  the  non-circular  mask(s)  used  in  conjunction  with  the 
convolution  error  term.  To  use  the  band-limit  mask  alone,  similar  to  the  runs  made 
in  experiment  two,  simply  answer  “1”  to  Q1  and  “0”  to  the  others. 

A. 4-3  Band-limit  Error  Term.  The  third  term  in  Jefferies  and  Christou’s 
error  function  is  a  band-limit  error  which  minimizes  the  high  frequency  content  of 
the  PSF.  This  error  term  was  used  in  experiment  three  with  mixed  results.  IDA 
adds  this  term  to  the  objective  function  whenever  the  telescope  aperture  parameter 
is  “2”.  Thus,  the  current  version  of  IDA  does  not  allow  for  use  of  the  band-limit 
error  without  employing  the  band-limit  mask  also.  As  explained  above,  entering  “2” 
for  Ql  will  require  additional  information  regarding  the  cutoff  frequency  of  each  PSF 
in  pixels.  The  band-limit  error  term  minimizes  all  pixel  values  in  the  OTF  estimate 
outside  of  1.39  times  the  input  cutoff  frequency  (12). 

As  explained  in  Section  4.4,  experiment  three  required  that  the  band-limit  error 
term  be  used  for  only  the  first  few  iterations,  then  IDA  was  restarted.  In  this  case, 
the  maximum  iterations  were  set  to  50  and  the  answers  to  Ql  through  Q4  were  2, 
0,  0,  0,  respectively.  After  running  50  iterations  with  the  band-limit  error,  IDA  was 
restarted  with  a  large  maximum  iteration  value  and  the  following  values  entered  for 
Ql  through  Q4:  1,  1,  0,  0.  The  first  “1”  makes  use  of  the  band-limit  mask  and  the 
second  “1”  tells  IDA  that  the  input  object  and  PSF  estimates  are  from  a  previous 
run  that  did  not  attain  a  minimal  objective  function  value. 

A. 4-4  Fourier  Modulus  Error.  The  third  question  (Q3)  allows  IDA  to 
employ  the  fourth  error  term  or  the  Fourier  Modulus  error.  This  error  term  was 
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not  explored  in  this  investigation  since  it  requires  explicit  knowledge  of  the  object 
attainable  only  through  other  image  processing  techniques.  The  Fourier  spectrum 
parameter  is  set  to  zero  unless  use  of  the  error  term  is  desired.  Entering  “2”  for  Q3 
requires  input  of  both  the  real  and  complex  portions  of  the  object  spectrum.  Simply 
entering  “1”  for  Q3  informs  IDA  that  Fourier  modulus  information  will  be  entered. 
Additionally,  use  of  the  Fourier  modulus  error  term  requires  the  input  of  a  “signal- 
to-noise”  filter  as  explained  in  Jefferies  and  Christou  [1993],  the  wavelength  of  the 
observations,  the  diameter  of  the  aperture,  and  image  scale  in  miUi-arcseconds/pixel. 
More  information  can  be  found  regarding  the  Fourier  modulus  error  in  Jefferies  and 
Christou  (12). 

A. 5  Modifications  and  Recompilation 

Since  the  code  supplied  by  Jefferies  and  Christou  is  in  FORTRAN,  array  size 
parameters  must  be  set  prior  to  compilation  of  the  code.  To  facihtate  processing 
of  images  having  different  sizes  or  support  areas,  a  header  file  containing  all  of  the 
applicable  parameters  for  processing  of  a  single  image  or  set  of  images  may  require 
modification  each  time  IDA  is  run.  This  header  file  sets  the  number  of  PSFs  that 
wiU  be  required  for  a  single  run,  the  number  of  pixels  on  one  side  of  the  image  array, 
the  number  of  variables  IDA  wiU  have  to  minimize  in  the  objective  function,  and 
the  tolerance  at  which  the  change  in  the  error  function  value  must  reach  before  a 
minimum  value  is  found.  Each  of  the  parameters  are  explained  in  detail  below. 

A. 5.1  Number  of  PSFs.  The  number  of  PSFs  or  npsf  must  be  set  to 
“1”  for  a  single  measurement,  object  and  PSF  combination.  However,  if  multiple 
realizations  of  the  same  object  are  input  into  IDA,  then  the  total  number  of  different 
PSF  estimates  which  IDA  will  have  to  make  at  each  iteration  is  the  value  for  npsf. 

A. 5. 2  Number  of  Pixels.  IDA  requires  that  all  arrays  input  into  the  pro¬ 
gram  (measurements,  estimates,  support,  etc.)  have  the  same  size,  npix  x  npix, 


A-6 


where  npix  represents  the  number  of  pixels  in  either  one  row  or  one  column  of  the 
square  array. 

A. 5. 3  Number  of  Variables.  This  parameter  gives  rise  to  modifying  the 
header  file  and  recompifing  IDA  even  for  different  images  of  the  same  array  size.  The 
number  of  variables,  nuar,  equals  the  total  number  of  non-zero  pixels  within  both  the 
object  and  PSF  support  arrays.  Therefore,  any  time  the  support  size  changes,  nvar 
changes  and  IDA  must  be  recompiled.  For  the  four  different  images  used  in  this 
investigation,  a  copy  of  the  header  file  was  kept  in  each  of  four  different  files  with 
the  appropriate  values  for  nvar.  When  a  run  was  required  for  image  B  for  instance, 
the  header  for  image  B  was  copied  into  the  header  file  using  the  following  UNIX 
command:  cp  header B  IDA.h.  Then,  compihng  was  accomphshed  through  a  “make¬ 
file”  supphed  by  Stuart  Jefferies  by  issuing  this  command:  make  IDA-FITS. 

A. 5. 4  Tolerance.  The  tolerance,  tol,  represents  a  very  small  amount  of 
change  in  the  value  of  the  error  function.  Once  the  amount  of  change  reaches  the 
tolerance,  then  it  is  assumed  that  the  conjugate  gradient  minimization  routine  has 
found  a  local  minimum  and  further  iterations  cease. 

A.  6  Summary 

This  appendix  provides  information  regarding  the  use  of  the  IDA  FORTRAN 
code  provided  by  Stuart  Jefferies  and  modified  for  use  in  this  investigation.  A  brief 
explanation  of  how  the  algorithm  is  implemented  was  presented  along  with  exphcit 
details  regarding  the  use  of  IDA  to  reconstruct  images.  The  program  allows  a  user 
to  utihze  many  different  means  of  implementing  the  unconstrained  minimization 
technique  to  solve  the  blind  deconvolution  problem.  It  is  hoped  that  this  appendix 
will  allow  wider  use  of  IDA  to  process  blurred  images  and  that  future  research  can 
proceed  more  efficiently  by  using  the  documentation  presented  here.  Copies  of  the 
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IDA  code  used  in  this  study  can  be  obtained  through  Major  Michael  Roggemann, 
Air  Force  Institute  of  Technology /ENP,  Wright- Patterson  AFB,  OH  45433-6583. 
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Appendix  B.  Simulated  Images 


B.  1  True  Objects  and  Noise-Free  Images 

This  section  contains  a  visual  representation  of  the  four  objects  used  in  exper¬ 
iments  one,  two  and  three.  Each  object  is  a  64  x  64  array  with  one  or  two  point 
sources  within  each  array.  Three  of  the  objects  contain  binary  stars  with  an  intensity 
ratio  of  1:0.7.  Immediately  after  each  “true”  object  appears  an  image  resulting  from 
the  convolution  of  the  object  with  a  single  PSF  shown  in  Figure  B.9. 
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Figure  B.9  The  true  PSF  was  used  to  generate  the  blurred  images  previously 
shown.  The  iterative  deconvolution  algorithm  attempts  to  estimate 
this  PSF  along  with  the  appHcable  object  for  each  blurred  image.  The 
blur  radius  for  the  PSF  shown  above  is  11  pixels  representing  the  spatial 
distance  from  the  center  of  the  array  to  the  first  zero. 


photo  events 


B.3  Photon  Limited  Images 


Image  A  (100,000  photo  events) 


Figure  B.IO  The  image  above  represents  a  short-exposure  of  Image  A  containing 
only  100,000  photo  events  {K  =  10®). 


Image  A  (10,000  photo  events) 


Figure  B.ll  The  image  above  represents  a  short-exposure  of  Image  A  containing 
only  10,000  photo  events  (A"  =  lO"*). 


Image  A  (5,000  photo  events) 


1'  a'.*  •'  *. 


Figure  B.12  The  image  above  represents  a  short-exposure  of  Image  A  containing 
only  5,000  photo  events  {^K  =  5  •  10^). 


Image  B  (100,000  photo  events) 


Figure  B.13  The  image  above  represents  a  short-exposure  of  Image  B  containing 
only  100,000  photo  events  {K  =  10®). 


Image  B  (10,000  photo  events) 


Figure  B.14  The  image  above  represents  a  short-exposure  of  Image  B  containing 
only  10,000  photo  events  {K  =  10^). 


Image  B  (5,000  photo  events) 


Figure  B.15  The  image  above  represents  a  short-exposure  of  Image  B  containing 
only  5,000  photo  events  {K  =  5  •  10®). 


Figure  B.16  The  image  above  represents  a  short-exposure  of  Image  C  containing 
only  100,000  photo  events  {K  =  10®). 


Image  C  (10,000  photo  events) 
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Figure  B.17  The  image  above  represents  a  short-exposure  of  Image  C  containing 
only  10,000  photo  events  (K  =  10'*). 


Image  C  (5,000  photo  events) 


Figure  B.18  The  image  above  represents  a  short-exposure  of  Image  C  containing 
only  5,000  photo  events  {K  =  5  •  10®). 
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Image  D  (100,000  photo  events) 


700 

600. 


Figure  B.19  The  image  above  represents  a  short-exposure  of  Image  D  containing 
only  100,000  photo  events  (K  =  10®). 


Image  D  (10,000  photo  events) 


80., 


Figure  B.20  The  image  above  represents  a  short-exposure  of  Image  D  containing 
only  10,000  photo  events  {K  —  10'*). 
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Image  D  (5,000  photo  events) 


Figure  B.21  The  image  above  represents  a  short-exposure  of  Image  D  containing 
only  5,000  photo  events  {K  =  5  •  10®). 
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Appendix  G.  Experiment  1  Results 

The  first  experiment  tests  the  unconstrained  minimization  approach  to  bhnd 
deconvolution  for  noise-free  blurred  images.  Results  are  presented  for  the  four  images 
shown  in  Figures  B.2,  B.4,  B,6  and  B.8  using  three  different  support  regions.  The 
support  regions  correspond  to  the  first,  second  and  third  zero  of  the  true  PSF  as 
explained  in  Section  3.7.2.  The  data  presented  in  this  appendix  refer  to  the  smallest 
support  size  (corresponding  to  the  first  zero)  as  support  region  1,  the  medium-sized 
support  as  support  region  2  and  the  largest  as  support  region  3.  A  ring  of  “spikes” 
appear  on  the  edge  of  the  support  region  for  each  object  estimate.  This  phenomenon 
seems  to  be  a  function  of  the  algorithm,  and  the  ampUtude  of  the  spikes  decreases 
with  larger  support  regions.  As  presented  in  Chapter  4,  the  best  results  are  obtained 
using  support  region  2  for  these  images. 
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Appendix  D.  Experiment  2  Results 


The  second  experiment  analyzes  the  performance  of  the  iterative  deconvolu¬ 
tion  algorithm  in  the  presence  of  photon  noise.  The  simulated  photon  Hmited  im¬ 
ages  shown  in  Figures  B.IO  through  B.21  are  used  as  input  measurements  into  the 
algorithm.  A  band-hmit  mask  in  the  convolution  error  is  compared  against  the  ba¬ 
sic  Lane  method  of  unconstrained  minimization.  Results  presented  in  Figures  D.l 
through  D.8  show  that  even  with  as  many  as  100,000  photo  events,  the  addition  of 
the  band-hmit  mask  greatly  enhances  the  “reconstructabihty”  of  the  algorithm.  Two 
PSF  estimates  are  shown  in  Figure  D.9  which  reveal  the  wide  variation  in  impulse 
response  estimation  in  the  presence  of  noise.  Presentation  of  further  PSF  estimates 
provide  httle  insight  to  the  abihty  of  this  algorithm  to  reconstruct  objects  from 
noisy  blurred  images,  therefore,  no  other  PSF  estimates  are  shown.  The  remainder 
of  the  data  presented  show  object  estimates  from  images  containing  10,000  and  5,000 
photo  events.  It  is  clear  that  less  light  reduces  the  algorithm’s  ability  to  effectively 
reconstruct  images. 


Object  A  Estimate 
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Figure  D.l  The  object  estimate  above  resulted  from  Image  A  [K  ~  10“)  using  the 
strict  Lane  method  for  unconstrained  minimization. 


Object  A  Estimate 
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Figure  D.2  The  object  estimate  above  resulted  from  Image  A  {K  =  10®)  using  a 
band-limit  mask  in  the  convolution  error  term. 


Object  B  Estimate 


Figure  D.3  The  object  estimate  above  resulted  from  Image  B  {K  —  10®)  using  the 
strict  Lane  method  for  unconstrained  minimization. 


Object  B  Estimate 
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Figure  D.4  The  object  estimate  above  resulted  from  Image  B  {K  =  10®)  using  a 
band-Hmit  mask  in  the  convolution  error  term. 


Figure  D.7  The  object  estimate  above  resulted  from  Image  D  {K  =  10®)  using  the 
strict  Lane  method  for  unconstrained  minimization. 


Object  D  Estimate 


Figure  D.8  The  object  estimate  above  resulted  from  Image  D  {K  =  10®)  using  a 
band-hmit  mask  in  the  convolution  error  term. 


PSF  Estimate 
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Figure  D.9  The  PSF  estimates  above  resulted  from  Images  A  and  D  (K  =  10®) 
using  a  band-limit  mask  in  the  convolution  error  term.  These  estimates 
reveal  that  the  presence  of  noise  in  the  measurement  causes  deviations 
in  the  PSF  estimate  from  the  “true”  PSF.  Subsequent  PSF  estimates 
are  not  shown. 


Object  A  Estimate 
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Figure  D.IO  The  object  estimate  above  resulted  from  Image  A  {K  =  10“*)  using  a 
band-limit  mask  in  the  convolution  error  term. 


Object  A  Estimate 


Figure  D.ll  The  object  estimate  above  resulted  from  Image  A  {K  =  5,000)  using 
a  band-hmit  mask  in  the  convolution  error  term. 
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Figure  D.12  The  object  estimate  above  resulted  from  Image  B  {K  =  10“*)  using  a 
band-limit  mask  in  the  convolution  error  term. 


Figure  D.13  The  object  estimate  above  resulted  from  Image  B  [K  =  5,000)  using 
a  band-Hmit  mask  in  the  convolution  error  term. 


Object  C  Estimate 


Figure  D.14  The  object  estimate  above  resulted  from  Image  C  {K  =  10^)  using  a 
band-limit  mask  in  the  convolution  error  term. 


Object  C  Estimate 


Figure  D.15  The  object  estimate  above  resulted  from  Image  C  {K  =  5,000)  using 
a  band-limit  mask  in  the  convolution  error  term. 
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Figure  D.16  The  object  estimate  above  resulted  from  Image  D  (if  =  10"*)  using  a 
band-bmit  mask  in  the  convolution  error  term. 


Figure  D.17  The  object  estimate  above  resulted  from  Image  D  (if  =  5,  000)  using 
a  band-bmit  mask  in  the  convolution  error  term. 


Appendix  E.  Experiment  3  Results 


The  third  experiment  tests  the  addition  of  a  band-limit  error  term  within  the 
objective  function  for  unconstrained  minimization  to  improve  the  reconstruction  of 
both  noise-free  and  photon  hmited  images.  The  band-limit  error,  as  explained  by  Jef¬ 
feries  and  Christou  (12),  is  utihzed  at  the  start  of  the  iterative  process  to  prevent  the 
trivial  solution  from  appearing.  Though  no  trivial  solution  was  previously  obtained 
with  these  images,  this  experiment  attempts  to  determine  whether  the  band-limit 
error  may  have  a  positive  effect  on  the  reconstruction  process.  This  experiment  used 
a  fixed  number  of  iterations  (50)  employing  the  band-limit  error  term,  then  ran  the 
algorithm  without  the  band-limit  error  term  to  achieve  proper  convergence.  Mixed 
results  appear  in  the  following  data  revealing  that  the  use  of  the  band-limit  term 
should  be  studied  further. 
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Figure  E.l  The  object  estimate  above  resulted  from  Image  A  (noise-free)  using 
band-Umit  error  term  in  the  unconstrained  minimization  technique. 


Figure  E.3  The  object  estimate  above  resulted  from  Image  A  {K  =  10“*)  using 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Object  B  Estimate 


Figure  E.5  The  object  estimate  above  resulted  from  Image  B  (if  =  10®)  using  a 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Object  B  Estimate 


Figure  E.6  The  object  estimate  above  resulted  from  Image  B  (K  =  10'^)  using  a 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Irradiance 


Figure  E.7  The  object  estimate  above  resulted  from  Image  C  (noise-free)  using 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Irradiance 


Figure  E.9  The  object  estimate  above  resulted  from  Image  C  (^  =  lO"*)  using 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Object  D  Estimate 


Figure  E.IO  The  object  estimate  above  resulted  from  Image  D  (noise-free)  using  a 
band-Umit  error  term  in  the  unconstrained  minimization  technique. 


Object  D  Estimate 


Figure  E.ll  The  object  estimate  above  resulted  from  Image  D  (if  =  10®)  using  a 
band-limit  error  term  in  the  unconstrained  minimization  technique. 


Object  D  Estimate 


Figure  E.12  The  object  estimate  above  resulted  from  Image  D  {K  =  10^)  using  a 
band-hmit  error  term  in  the  unconstrained  minimization  technique. 
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